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ABSTRACT 


Liquefaction  causes  a large  portion  of  all  damage  done  by  earthquakes.  The  damage 
is  especially  severe  to  lifeline  structures  such  as  pipelines.  This  report  examines  the  state-of- 
the-art  of  the  application  of  System  Identification  (SI)  methods  to  the  liquefaction  problem, 
with  special  attention  to  lifelines.  System  identification  is  seen  as  the  best  way  to  ascertain 
large  strain  soil  properties  in  situ.  A thorough  introduction  to  SI  methods  and  spectral 
analysis  is  given.  The  traditional  Fourier-based  methods  are  found  to  be  inexact  since  the 
sample  variance  is  equal  to  the  sample  mean  if  averaging  techniques  are  not  used.  There 
is  an  additional  problem  since  earthquake  signals  are  not  stationary.  Autoregressive-moving 
average  models  are  seen  as  a better  analysis  method,  especially  the  newer  adaptive  methods 
that  are  designed  for  non-stationary  signals.  A significant  bibliography  is  included. 
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CHAPTER  1 INTRODUCTION 


1.1  Background 

1.1.1  Earthquakes  and  Liquefaction 

Over  the  years,  some  of  the  most  spectacular,  and  costly  damage  caused  by 
earthquakes  has  been  due  to  liquefaction  of  sands.  A functional  description  of  liquefaction 
can  be  given  quite  simply.  When  a suitably  intense  earthquake  shakes  a loose,  saturated 
sand,  the  grains  tend  to  consolidate  into  a more  compact  packing.  Since  all  these 
movements  are  happening  very  rapidly,  there  is  no  chance  for  the  volume  to  reduce  through 
pore  water  dissipation.  The  incompressible  pore  fluid  then  takes  up  all  the  applied  stress, 
the  effective  stress  goes  to  zero,  and  the  deposit  "liquefies."  Since  a liquid  has  no  shear 
strength,  disastrous  consequences  occur. 

The  disastrous  consequences  of  liquefaction  was  brought  to  the  fore  in  1964  by  the 
Niigata  and  Alaska  earthquakes  of  that  year.  Liquefaction  also  triggers  earth  slides  and 
large  displacements  of  earthen  dams.  A large  slide  caused  by  the  1970  Peruvian  earthquake 
killed  over  18,000  people  (EERI,  1986).  A terrible  disaster  was  narrowly  avoided  when  the 
San  Fernando  dam  suffered  very  large  displacements  due  to  the  1971  San  Fernando 
earthquake.  Similar  damage  has  occurred  over  the  years  in  locations  as  diverse  as  China, 
Nicaragua,  Japan,  Charleston,  SC,  San  Francisco,  the  Imperial  Valley,  CA,  and  Idaho. 

The  effects  of  liquefaction-caused  damage  to  lifelines  are  especially  costly.  Damage 
to  roads,  rail,  telecommunications,  power,  and  pipelines  of  all  types  is  always  harmful,  but 
is  especially  so  during  time  of  emergency.  One  of  the  most  striking  examples  of  the  effect 
of  lifeline  damage  on  public  safety  is  the  occurrences  during  the  1906  San  Francisco 
earthquake.  After  that  earthquake,  over  490  city  blocks  were  totally  destroyed  by  a fire,  the 
largest,  most  deadly  fire  in  U.S.  history  (O’Rourke  et  al.,  1991).  Little  could  be  done  to  stop 
the  spread  of  the  fire  since  the  pipelines  carrying  water  were  broken  due  to  liquefaction- 
induced  ground  displacements.  It  was  estimated  that  56  percent  of  the  municipal  water 
supply  was  completely  cut  off. 


1.1.2  Estimation  Techniques  and  Liquefaction 

The  reader  might  wonder,  what  does  system  identification  and  estimation  of  large 
strain  soil  properties  have  to  do  with  liquefaction?  While  not  immediately  obvious,  there 
are  several  important  reasons  to  pay  attention  to  these  methods. 

The  most  important  piece  of  knowledge  to  be  gained  is  that  concerning  the  actual 
behavior  of  soils  under  strong  motion  excitation.  At  present,  the  modulus  degradation  and 
effective  damping  ratio  curves  are  based  only  on  laboratory  tests.  Laboratory  tests  will  all 
yield  at  best  approximate  results  since  no  one  can  run  a laboratory  test  on  an  undisturbed 


loose  sand.  The  preliminary  reports  from  back  calculating  earthquake  response  imply  that 
the  laboratory  degradation  curve  might  be  too  high  at  intermediate  strains.  In  addition,  the 
results  from  two  independent  methods  (Chang  et  al.,  1990;  Abdel-Ghaffar  and  Scott,  1979) 
show  that  the  customary  hyperbolic  shape  of  the  laboratory  damping  curve  might  be 
incorrect,  and  actually  is  S-shaped. 

Knowledge  about  the  interplay  between  pore  water  pressure  build-up,  strength  of 
shaking,  and  soil  nonlinearity  is  fundamental  to  rational  design  for  earthquake  loading.  This 
interplay  can  only  be  studied  by  an  analysis  of  undisturbed  real-life  situations  using  complete 
sets  of  data  such  as  available  from  Lotung.  It  is  of  utmost  importance  for  the  database  to 
be  enlarged  with  results  from  other  well-instrumented  site.  Unfortunately,  there  is  a grave 
lack  of  such  sites. 

The  response  of  soil  to  strong  motion  can  be  used  for  site  characterization  and 
microzonation  (Finn,  1991).  In  this  case,  the  amplification  factor,  or  spectral  ratio,  is  the 
important  parameter  being  sought.  The  amplification  factor  is  just  the  transfer  function, 
which  can  be  estimated  much  more  accurately  using  the  system  identification  methods.  A 
simplified  microzonation  analysis  of  Charleston,  South  Carolina  was  made  using  SHAKE  as 
the  analysis  tool  (Elton  and  Martin,  1989),  and  the  amplification  effects  of  geological 
structures  have  been  examined  theoretically  (e.g.  Faccioli,  1991)  and  experimentally  (e.g. 
Silva,  1989;  Bard  and  Gabriel,  1986). 

A final  incentive  to  study  in  situ  behavior  of  soils  during  earthquake  strong  motion 
is  brought  up  in  a late  paper  by  Rolhns  and  Seed  (1990).  This  is  the  question  of  what 
influence  structures  have  on  their  founding  soils.  In  particular,  does  a structure  increase  or 
decrease  the  liquefaction  potential,  and  if  the  soil  does  liquefy,  will  the  resulting 
displacement  be  more  or  less  than  the  free-held?  These  are  important  questions  since  the 
present  method  of  analysis  and  design  implies  a free-held  condition.  The  available  data  is 
quite  scanty,  based  mostly  on  a few  shaking  table  and  centrifuge  tests. 

Rollins  offers  some  conclusion  as  to  the  effects  of  different  types  of  structures.  These 
include  indications  that  excess  pore  pressure  ratios  might  be  signihcantly  lower  beneath  a 
structure,  and  that  the  soil  near  a foundation  is  more  susceptible  to  generation  of  excess 
pore  pressure  than  the  free-held.  Free  held  analysis  appears  to  be  too  conservative  for  long- 
period  structures  on  medium-dense  sands,  and  unconservative  for  short-period  buildings  on 
loose  to  medium  sands. 

In  a wide-ranging  paper,  Ambraseys  and  Sarma  (1969)  give  calculations  showing  that 
concentrated  loads  from  a structure  "may  cause  much  more  widespread  liquefaction  effects 
than  local  inhomogeneities."  This  is  through  local  failure  acting  as  seed  for  widespread 
progressive  failure.  On  the  other  hand,  the  local  increased  vertical  effective  stress  due  to 
the  structure  will  decrease  the  chance  of  liquefaction.  Finding  the  balance  between  the  two 
effects  leads  to  the  same  questions  as  Rollins  and  Seed. 
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The  ideal  way  to  find  definitive  conclusions  to  these  important  questions  is  to  utilize 
input  and  output  soil  motions  recorded  in  the  free-held  and  adjacent  to  a structure.  This 
data  is  available  right  now  for  one  case  — the  Lotung  site.  A detailed  analysis  using  state-of- 
the-art  system  identification  techniques  is  needed. 


1.1.3  Purpose 

This  report  was  written  to  evaluate  the  current  state-of-the-art  of  in  situ  methods  of 
soil  property  measurement,  which  allow  accurate  prediction  of  liquefaction  potential,  and 
the  possible  displacements  if  liquefaction  did  occur.  The  report  summarizes  and  evaluates 
significant  technical  papers  on  the  use  of  system  identification  methods  for  estimating  soil 
parameters  needed  to  understand  large  strain  soil  behavior  during  earthquake  loading.  This 
topic  is  of  direct  import  to  the  behavior  of  lifelines.  This  report  makes  no  attempt  to 
identify  and  enumerate  every  paper  or  technical  publication  written  on  these  subjects.  It 
serves,  rather,  as  a thorough  overview  and  evaluation  of  where  the  profession  is  today. 

The  report  assumes  some  degree  of  technical  sophistication  by  the  reader,  although 
an  attempt  is  made  to  explain  complicated  or  unfamiliar  material.  The  liberal  use  of 
references  should  allow  the  reader  to  find  an  understandable  source  of  explanation  for  most 
topics  discussed.  The  report  often  takes  a "critical"  point  of  view  when  examining  proposed 
methods.  This  is  to  examine  the  underlying  suppositions  made  by  a given  approach,  which 
define  the  validity  and  applicability  of  that  method. 


1.2  Scope 

This  report  examines  the  field  of  system  identification,  and  its  applicability  to 
estimating  the  behavior  of  soils  undergoing  strong  motion.  Since  the  field  of  system 
identification  is  unfamiliar  to  most  geotechnical  engineers.  Chapter  2 gives  an  introduction 
to  the  meaning  of  relevant  techniques.  Of  great  importance  is  the  thorough  evaluation  of 
the  implications  of  spectral  analysis. 

Chapter  3 examines  pertinent  applications  of  system  identification  to  liquefaction-type 
problems.  Special  attention  is  given  to  using  system  identification  methods  to  estimate 
mechanical  properties  of  soil  undergoing  intermediate-to-large  strains,  which  can  not  be 
examined  by  geophysical  in  situ  techniques. 

Chapter  4 contains  the  important  conclusions  reached  in  this  report,  and  briefly 
summarizes  key  concepts.  A very  thorough  bibliography  is  documented  in  the  final  chapter. 

Finally,  it  should  be  mentioned  that  the  most  thorough  and  complete  work  on  the 
subject  of  liquefaction  is  the  report  written  by  the  Committee  on  Earthquake  Engineering 
of  the  National  Research  Council,  1985. 
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CHAPTER  2 PRINCIPLES  OF  PROCESS  CHARACTERIZATION  IN 
THE  TIME  AND  FREQUENCY  DOMAINS 


2.1  Estimates  of  Soil  Properties  from  Dynamic  Behavior 

2.1.1  Introduction 

A major  component  of  the  lifeline  infrastructure  are  buried  pipelines.  Over  the  years, 
analytical  computer  programs  have  been  developed  to  calculate  the  effects  of  earthquakes 
on  these  lifelines.  As  the  programs  are  able  to  make  more  and  more  accurate  predictions, 
the  significance  of  the  accuracy  and  relevancy  of  the  input  data  becomes  more  and  more 
important.  Most  of  these  analyses  use  a Winkler  model,  which  replace  the  soil  with  springs 
and  dashpots  (Zhang  et  al.,  1991).  The  relevant  soil  parameters  become  dynamic  stiffness 
and  damping.  Since  earthquakes  excite  the  soil  past  the  elastic  range,  estimates  of  the 
strain-dependant  soil  stiffness  and  damping  are  essential  for  further  improvements  in  lifeline 
design. 


Field  geophysical  techniques  are  relevant  to  liquefaction  analysis  in  so  far  as  they 
provide  the  shear-wave  velocities  of  a site.  The  velocities  can  easily  be  converted  to  soil 
stiffness,  or  correlations  with  liquefaction  potential  derived  directly  for  S-wave  velocity. 
However,  the  shear  modulus  thus  calculated  is  the  small  strain  modulus,  or  and  is  only 
valid  for  the  elastic  region  of  the  soil.  The  limitation  is  due  to  the  inability  to  reliably  impart 
strains  into  the  soil  much  greater  than  IxlO"^.  Therefore,  it  has  been  impossible  to  measure 
threshold  strain,  y-p,  and  the  modulus  degradation  curve,  G/Gjjjax^  in  situ.  Direct  knowledge 
of  these  variables  would  bring  the  ability  to  make  quantitative  estimates  of  potential  and 
possible  displacements  closer  to  reality.  Another  important  benefit  would  be  the  possibility 
to  measure  nonlinear  material  damping. 

Attempts  to  input  enough  energy  into  the  ground  to  cause  intermediate  to  large 
strains  have  not  been  very  successful.  The  amount  of  energy  needed  would  destroy  a bore- 
hole, and  would  be  destructive  on  the  surface  as  well.  There  is  also  the  problem  of  the 
transducers  being  in  the  near-field  if  they  are  close  to  a source  large  enough  to  cause  large 
strains  in  an  immediate  area.  One  exception  was  a project  undertaken  for  the  Nuclear 
Regulatory  Commission  (Shannon-Wilson,  1976)  where  intermediate-to-large  strains  were 
input  in  a large  scale  field  experiment. 

An  obvious  example  of  large  strain  experiment  would  be  the  use  of  high  explosives. 
While  explosives  have  been  used  in  the  Soviet  Union  to  estimate  liquefaction  parameters 
(Florin  and  Ivanov,  1961),  it  was  done  to  develop  a correlation  with  settlement,  and  no 
geophysical  measurements  were  made.  Positive  results  of  a large  in  situ  impulse  test  were 
reported  (Shannon-Wilson,  1976).  However,  there  has  been  no  follow-up  on  this  work,  and 
other  researchers  have  not  attempted  similar  studies.  It  is  not  known  at  the  present  time 
whether  that  is  because  of  lack  of  interest  or  problems  with  the  reported  method. 


The  optimum  situation  would  be  the  ability  to  make  measurements  during  different 
magnitudes  of  earthquake  excitement.  In  this  case  shear  strain  in  the  layers  of  interest,  and 
stiffness  (velocity),  would  be  continually  monitored.  Since  earthquakes  can  not  be  made-to- 
order,  the  chances  of  this  situation  happening  are  virtually  nonexistent.  The  instrumentation 
would  also  be  extremely  difficult.  However,  use  of  inverse  theory  allows  the  soil  parameters 
of  interest  to  be  calculated  from  attainable  data  — the  ground  motion  records  of  the  motion 
going  into  the  layers  of  interest,  and  above  the  layer  itself.  This  set  of  information  is 
available  for  several  sites  (Chang  et  al.,  1991,  1990;  Holzer  et  al.,  1989). 

The  modeling  of  a mechanical  system  as  a transfer  function  calculated  from  known 
a input-output  history  is  commonly  called  system  identification  (SI).  If  a suitable  model  is 
chosen  to  represent  the  system  of  interest,  the  model  parameters  derived  will  correspond 
to  important  mechanical  parameters  of  the  system,  such  as  damping,  natural  frequency,  and 
stiffness.  Often  times,  SI  is  the  only  method  available  to  estimate  these  properties,  especially 
since  the  method  does  not  actively  interfere  with  the  material  properties  being  measured. 
This  chapter  will  investigate  the  theory  required  to  undertake  this  system  identification  for 
liquefiable  soils. 


2.1.2  Modeling  a System 

In  the  simplest  case,  a layer  of  soil  can  be  modelled  as  a linear  system,  as  diagramed 
in  Fig.  2.1.  In  the  time  domain,  the  filtering  process  of  a signal  passing  through  this  layer 
is  represented  as  convolution,  Eq.  2.1: 

n-l 

yin)=Yl  x(myh(n-m)=x*h,  (2-1) 

OT=0 


where  y, 

Xt 

ht 


n 

♦ 


= output  time  series, 

= input  time  series, 

= filter  vector  or  impulse  response  function, 
= 1... length  of  y(t), 

= convolution. 


The  process  of  inversion,  or  deconvolution,  allows  the  estimation  of  the  system  response 
function  (filter)  if  the  input  and  output  signals  are  known.  Theoretically,  the  input  and 
output  vectors  represent  the  coefficients  of  a polynomial  (Z-transform)  and  the  system 
response  function  h^  can  be  solved  for  by  polynomial  division  (Bracewell,  1978).  However, 
if  there  is  any  noise  present,  and  there  ALWAYS  is  (due  to  quantization  error  if  nothing 
else),  the  quotient  is  irrational  and  frequently  becomes  unstable. 

The  usual  method  of  time  domain  deconvolution  is  the  least  squares  approach  (Silvia 
and  Robinson,  1979).  While  the  deconvolution  process  is  non-unique,  the  least  squares 
method  yields  a system  filter  that  is  unique  in  a mathematical  (least  squared  error)  and 
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physical  sense  (minimum  phase).  There  are  other,  more  involved,  time  domain  solutions 
which  are  not  in  general  use  (e.g.  Simmons,  1991).  The  time  domain  schemes,  while 
possessing  certain  advantages,  are  very  computation  intensive.  The  common  solution  is  to 
deconvolve  in  the  frequency  domain. 

In  the  frequency  domain,  convolution  becomes  a simple  multiplication  (Bracewell, 
1978),  and  Eq.  2.1  transforms  into  Eq.  2.2: 


(2.2) 


y(cD)=X(a))*/f(co), 


where  o)  = circular  frequency. 


= frequency  function  of  input  time  series, 

= frequency  function  of  output  time  series, 

Hy  = system  frequency  response  function. 

For  this  report,  time  domain  functions  will  be  indicated  by  lower  case  variables  and 
frequency  domain  functions  by  upper  case  variables,  with  t and  (■>  dropped  when  obvious. 
The  transformation  into  the  frequency  domain  speeds  computation  but  does  not  diminish 
the  impact  of  noise  on  the  calculation.  If  the  variance  of  the  measurement  is  close  to  the 
same  order  of  magnitude  as  the  signal  energy  at  a given  frequency,  the  results  of  the  simple 
calculation  shown  in  Eq.  2.2  are  seriously  flawed  (Newland,  1984). 

2.2  A Traditional  Approach  to  System  Identification 
2.2.1  Introduction 

As  in  the  time  domain,  there  is  an  acceptable  frequency  domain  alternative  method 
for  computing  the  frequency  response  function  using  the  auto  and  cross  spectrum.  Using 
these  functions,  the  frequency  response  function  is  defined  strictly  for  stationary  input  signals 
as  (Bendat  and  Piersol,  1986): 


(2.3) 


G,  (XX-) 


where  = auto-spectrum  = XX*  = Fourier  transform  of  R^x 

Gxy  = cross-spectrum  = YX*  = Fourier  transform  of  R^y 

Rx^  = autocorrelation  of  x 
Rxv  = cross-correlation  of  x with  y 


(2.3a) 

(2.3b) 


= complex  conjugate. 
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The  use  of  the  cross-spectrum  in  calculating  the  frequency  response  function  is  little  more 
complicated  than  Eq.  2.2,  but  has  the  added  benefit  that  the  statistical  variability  of  is 
canceled  by  the  variability  of  The  real  part  of  H is  called  the  system  gain  factor,  while 
the  phase  information  is  carried  by  the  imaginary  part.  The  frequency  response  of  the 
system  is  now  experimentally  characterized  by  H,  but  useful  mechanical  parameters  such  as 
damping  and  sti^ess  are  still  unknown.  The  key  is  to  link  the  frequency  response  function 
to  the  defining  equations  of  the  system. 


2.2.2  Modeling  Simple  Mechanical  Systems  - Force  Input  and  Displacement  Output 

The  simple  layered  soil  system  in  question  can  be  characterized  as  a damped  single- 
degree-of-freedom  (SDOF)  system.  The  most  basic  dynamic  situation  is  with  a known  input 
force,  and  the  resultant  displacement  measured,  as  pictured  in  Fig.  2.2.  While  this  degree 
of  knowledge  is  unrealistic  for  the  earthquake  system  in  question,  it  is  a good  model  for 
surface  excitation  methods.  The  forces  acting  on  the  mass  are  given  by  Eq.  2.4: 

(2-4) 


and  the  equation  of  motion  is 


nty+cy+ky=f. 


(2.5) 


where  f(t) 

= input  force 

fk(t) 

= -ky 

= spring  force 

fc(t) 

= -cy 

= damping  force) 

f.(t) 

= -my 

= inertial  force 

y 

= output  system  displacement 

y 

= dy/dt 

= velocity 

y 

= dVdt^ 

= acceleration. 

The  system  response  function  is  defined  for  an  impulsive  forcing  function.  The  Fourier 

transform  of  the  impulse  response  function,  Y,  is  the  frequency  response  function,  H.  Based 
on  the  definition  of  the  Fourier  transform  and  related  theorems  (Bracewell,  1978)  the 
transform  of  Eq.  2.5,  with  f(t)  = becomes 

(-Ci)^m+i(oc+it)T=l 


where  i =7-1 

and  the  frequency  response  function  is 

Equation  2.7  is  simplified  by  defining  the  damping  ratio 
and  the  undamped  circular  natural  frequency 


8 


(2.7) 


H=Y=- 


1 


c=- 


2^/i^ 


(2.8) 


(2.9) 


Equation  2.7  is  now  written  as 


H = 


1 

k 


1- 


( 

+ i2C 

) 


/ 

h) 


(2,10) 


For  ease  of  analysis,  the  frequency  response  function  can  be  broken  into  a system 
gain  factor  and  a phase  factor  by  writing  Eq.  2.10  in  polar  form.  The  gain  factor  is  defined 
as 


\ 

k 


\ 


1- 


\fn, 


12 


2C 


\fn/ 


(2.11) 


and  the  phase  factor 


(J)=taii  ^ 


2C 


1- 


-T 

\-^«j  j 


(2.12) 


The  gain  factor  is  shown  graphically  in  Fig.  2.3.  The  gain  factor  for  this  particular  forced 
system  is  called  the  magnification  function.  Note  that  for  zero  frequency  the  magnification 
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function  is  the  inverse  of  the  spring  constant  k.  The  resonant  frequency  fj.  can  be  read  from 
the  graph  of  the  gain  factor,  while  the  fundamental  frequency  4 always  corresponds  to  a 
phase  factor  of  1.57  The  damping  ratio  can  be  calculated  from  Eq.  2.13: 


Sr 

/. 


2 


(2.13) 


or  by  a simple  half-power  graphical  technique  (Richart  et  al.,  1970).  All  the  system 
parameters  can  be  derived  from  the  plots  of  the  frequency  response  function  H calculated 
by  Eq.  2.3. 

2.2.3  Modeling  Earthquake  Response  — Displacement  Input  and  Displacement  Output 

The  field  problem  being  addressed  is  that  of  a soil  layer  where  the  displacement  into 
the  layer  (the  bottom)  and  the  displacement  at  the  top  of  the  layer  are  known.  This  system 
is  modeled  by  the  mechanical  analog  shown  in  Fig.  2.4.  The  goal  is  to  use  system  parameter 
identification  techniques  to  compute  the  soil  stiffness  and  damping  from  input  and  output 
records,  as  in  the  previous  example. 

The  equation  of  motion  for  this  system  is: 

my^cy+ky  = kx+cx.  (2-14) 


All  the  variables  are  defined  as  in  the  previous  example.  Applying  the  Fourier  transform 
to  Eq.  2.14  yields  the  frequency  domain  equation  Eq.  2.15: 

(-a)^/w+i(oc+A:)y'=(Jt+zci)c)  (2-15) 

and  the  frequency  response  function  H is  defined  as 

— . (2.16) 

k-isp-m+imc 


10 


Xt  ^ 

Constant  Parameter  Linear  System 

Yt 

X, 

h, 

Ho, 

Yco 

Fig.  2.1  Diagram  of  a linear  system  (filter). 


Fig.  2.2  Single-Degree-of-Freedom  system,  force  input-displacement  output  (Bendat  and  Piersol,  1986). 


Fig.  2.3  Frequency  response  function  (Gain  Factor)  of  SDOF  system  with  force  input  (Bendat  and  Piersol,  1986). 
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Substituting  damping  ratio  (Eq.  2.8)  and  fundamental  frequenq^  (Eq.  2.9)  into  Eq.  2.16  yields 


l+i2C 


£ 

\fn, 


1- 


+ i2C 


'£ 

/rtj 


(2.17) 


The  gain  factor  of  H,  defined  in  Eq.  2.18,  is  called  the  transmissibility  function  since 
it  defines  how  much  of  the  forcing  displacement  is  transmitted  through  the  system. 


1+ 

2C 

/) 

/J 

2 

( fV 
1-  ^ 

12 

+ 

2C 

f/)? 
£/ . 

(2.18) 


The  phase  factor  is  given  by: 


4)  =taii  1 


2C 


1- 


f/ 

\fnf 


\2 


+ 4C^ 


v-^«/ 


(2.19) 


The  transmissibility  function  is  shown  in  Fig.  2.5. 

While  the  graph  of  the  gain  factor  for  the  displacement-displacement  system  look 
similar  to  that  of  the  force-displacement  system  (Fig.  2.3),  the  equations  are  very  different. 
From  Fig.  2.5  it  is  seen  that  for  any  system  parameters,  the  gain  factor  equals  one  at  a 
frequency 

(i)  = v^*(»)^  (2.20) 

The  maximum  value  of  the  gain  factor  is  related  to  the  damping  ratio  by 

|g|=  ^ (2.21) 

\/l6C‘'-8C’-2+2yi+8C^ 

which  is  shown  in  graphical  format  in  Fig.  2.6  (Crede,  1957).  For  this  system,  values  can 
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Fig.  2.4  Single-Degree-of-Freedom  system,  displacement  input-displacement  output 
(Bendat  and  Piersol,  1986). 


Fig.  2.5  Frequency  response  function  of  SDOF  system  with  foundation  motion  input  (Crede,  1951). 
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Fig.  2.6  Graphical  relationship  between  Transmissibility  and  dannping  ratio  for  SDOF  system 
with  foundation  motion  input  (Crede,  1951). 
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only  be  derived  for  the  damping  ratio  and  resonant  frequency.  There  is  not  enough 
information  to  directly  identify  the  actual  system  stiffness  and  damping.  In  this  case,  mass 
is  an  unknown  that  prevents  identification  of  these  parameters.  The  concept  of  parameter 
mass  is  not  straight-forward  for  a layer  of  soil.  However,  it  is  quite  possible  that  the  use  of 
the  mass  of  a unit  area  column  of  soil  will  provide  a correct  solution. 


2.2.4  Evaluation  of  Simple  System  Identification  Techniques 

The  inability  to  define  all  the  system  parameters  highlights  an  important  difficulty  of 
the  system  identification  method.  Even  for  a simple  SDOF  system  with  realistic  boundary 
conditions  there  is  not  a closed-form  solution  for  all  the  basic  parameters.  For  more 
complicated  systems  the  equations  get  even  more  intractable  and  arcane.  Often  times  it  is 
not  possible  to  find  a suitable  equation  to  define  the  real-world  system  of  interest. 

The  solution  given  above  is  not  even  valid  for  the  real-life  situation  of  earthquake 
loading  of  a soil  profile.  The  derivation  is  strictly  valid  for  a time  invariant  system,  which 
does  not  accurately  describe  a soil  strained  past  threshold  strain.  The  stiffness  and  damping 
of  the  soil  are  not  linear,  and  measuring  this  non-linearity  is  one  goal  of  utilizing  system 
identification.  In  addition,  the  earthquake  forcing  function  is  not  time-stationary,  so  common 
spectral  estimation  methods  cannot  be  used  to  reduce  effects  of  noise  without  extensive 
increase  in  the  difficulty  of  the  computation.  The  non-stationary  methods  solve  for  the 
frequency  response  function  for  input  excitations  that  vary  through  time  (Bendat  and  Piersol, 
1986;  Newland,  1984).  For  the  traditional  "Bendat  and  Piersol"  approach,  this  involves  using 
energy  spectral  density  estimates  rather  than  power  spectral  density  estimates  since  the 
period  T is  finite  rather  than  infinite  time  for  which  the  power  spectrum  is  defined.  The 
single  summations  of  Eqns.  2.3  become  double  summations,  since  the  effect  of  change 
through  time  must  be  actively  accounted  for.  The  double  sum  estimates  are  not  as  robust 
as  the  ones  for  time  invariant  systems  since  the  estimation  are  made  for  two  stochastic 
variables  with  cumulative  variances. 

Assuming  for  the  moment  a time  invariant  system,  the  question  is  whether  the 
excitation  is  stationary  or  not.  If  the  input  is  stationary,  then  the  cross-spectral  technique 
can  be  used  to  define  the  frequency  response  function,  minimizing  the  effects  of  noise. 
Stationary  excitation  is  only  to  be  expected  for  some  ambient  vibration  tests  or  from 
controlled  oscillatory  sources.  For  these  two  cases  the  strain  levels  are  small  and  the  soil  will 
behave  linearly.  The  forcing  function  of  the  oscillatory  source  is  known  so  that  all  system 
parameters  can  be  estimated  using  the  force  displacement  model.  If  the  masses  of  the 
system  are  known,  the  problem  can  be  solved  uniquely  for  a multiple  layered  soil  profile 
(Udwadia,  1985). 

Earthquake  excitation  is  obviously  non-stationary.  It  is  an  unique  event  changing 
through  time  that  can  not  ever  be  repeated  so  the  effects  of  noise  can  not  be  eliminated 
through  direct  averaging.  Because  the  input  is  non-stationary  the  power  spectrum  methods 
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are  strictly  invalid  since  they  depend  on  a statistically  representative  sampling  of  an  infinite 
length  signal.  However,  claim  has  been  made  that  the  spectral  methods  are  equally  valid 
for  non-stationary  transient  data  if  there  are  a sufficient  number  of  trials  so  that  a valid 
expected  value  of  the  auto-  and  cross-correlations  can  be  calculated  (Bendat,  1990;  Bendat 
and  Piersol,  1986).  This  ensemble  averaging  serves  to  improve  the  statistical  certainty  of  the 
estimates. 


2.3  Estimates  of  a Process  in  the  Frequency  Domain 
2.3.1  Fourier-based  Methods 

At  this  point  it  is  important  to  examine  what  exactly  is  meant  by  the  spectral 
representation  of  a time  series,  and  the  results  of  the  various  methods  of  calculating  the 
"spectrum."  The  most  familiar  form  of  spectral  estimation  is  based  on  the  Fourier  transform 
of  a continuous  signal. 


F((o)=ffit)€-'^‘dt.  (2-22) 


This  estimate  is  non-parametric  since  no  particular  model  of  process  is  assumed  in  the 
formulation  of  the  estimate.  However,  the  process  is  arbitrarily  assumed  to  be  representable 
by  an  orthogonal  basis  function  of  harmonic  sinusoids.  The  basis  function  could  just  as  easily 
be  assumed  to  be  a linear  combination  of  damped  exponentials  (Prony’s  method),  a finite 
number  of  arbitrary  complex  sinusoid  in  white  noise  (Pisarenko  harmonic  decomposition), 
or  the  output  of  a sharp  bandpass  filter  centered  at  each  ft'equency  of  interest  (Capon’s 
Maximum  Likelihood  Method)  (see  Kay  and  Marple,  1981). 

Since  F(o)  is  a complex  value,  the  following  relation  exists: 

F(a))=|F(a))|  (2-23) 

where  | F(ci))  | = amplitude  spectrum 

gi<^)(«)  _ pj^ase  delay  spectrum. 

The  conservation  of  energy  between  the  time  and  Fourier  domain  is  given  by  Parseval’s 
energy  theorem  in  terms  of  the  square  of  the  amplitude  spectrum  (Bracewell,  1978).  The 
energy  spectral  density  is  defined  as  the  square  of  the  amplitude  spectrum,  | F(ci>)  | and 
gives  the  distribution  of  energy  as  a function  of  a set  of  harmonic  sinusoid.  The  square  of 
the  amplitude  spectrum  is  commonly  called  the  periodogram  or  power  spectrum  (Robinson, 
1982).  It  is  related  to  the  autocorrelation  of  the  same  time  function  by  the  Wiener- 
Khintchine  theorem  which  states  that  the  power  spectrum  is  the  Fourier  spectrum  of  the 
autocorrelation. 
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In  practice,  the  data  is  available  for  a finite  time  period,  and  the  data  is  in  discrete 
form  for  manipulation  by  a digital  computer.  This  necessitates  the  development  of  the 
discrete  Fourier  transform  (DFT)  (Brigham,  1974): 


F 


f— ] 

UA(j 


N-l 


=Y:xikAt) 

*=0 


(2.24) 


Where 


n 

N 

T 

Af 

At 


= 0,  1,  2,  ...,  N-l  = discrete  frequency  counter 

= number  of  samples 
= digitized  period  = NAt 

= frequency  resolution  1/NAt 
= time  between  samples. 


The  amplitude  spectrum  and  the  energy  spectral  density  are  defined  as  for  the  continuous 
case.  However,  the  values  for  the  discrete  case,  taken  at  the  preassigned  regular  Fourier 
frequencies,  do  not  match  those  for  the  analog  function.  This  is  because  the  sampled 
spectrum  is  actually  the  product  of  F(ci>)  with  the  sine  function  (transform  of  the  time  domain 
boxcar  windowing  function  which  is  convolved  with  the  signal  during  sampling).  "Thus  the 
discrete  spectrum  based  on  a finite  data  set  is  a distorted  version  of  the  continuous  spectrum 
based  on  an  infinite  data  set."  (Kay  and  Marple,  1981). 


The  unavoidable  windowing  of  the  data  in  the  time  domain,  and  multiplication  of  the 
spectrum  with  the  transform  of  the  window,  leads  to  a smearing  of  the  data  referred  to  as 
leakage.  For  the  simplest  example  of  the  boxcar  window,  the  DFT  representation  of  a 
simple  sinusoid  will  be  spread  out  over  a main  lobe  width  proportional  to  1/NAt  with  ripples 
occurring  at  discrete  intervals.  These  ripples  are  a function  of  the  mathematics  alone,  and 
are  not  present  in  the  actual  infinite  length  analog  signal.  Use  of  other  window  shapes  can 
decrease  the  amount  of  ripple,  but  at  the  expense  of  widening  the  main-lobe,  thus  decreasing 
the  frequency  resolution  of  the  transform.  The  various  windows  used  also  have  the 
unwanted  effect  of  biasing  the  data  towards  the  time  associated  with  the  peak  of  the  window 
(Geckinli  and  Yavuz,  1978). 


The  windowing  of  data  is  based  on  the  assumption  that  the  portions  of  the  times 
series  outside  the  window  are  zero.  For  the  case  of  a transient  such  as  an  earthquake 
seismogram,  this  assumption  is  realistic  and  bias  is  held  to  a minimum.  For  an  ongoing 
process  such  as  due  to  a mechanical  oscillator,  the  assumption  is  unrealistic  and  the  resulting 
spectrum  will  be  smeared  and  biased.  For  narrow-band  signals,  the  overlapping  of  adjacent 
side-lobes  can  hide  the  existence  of  close-by  lower-energy  components.  However,  for  wide- 
band processes  leakage  is  not  such  an  important  problem.  For  a soil-structure  interaction 
problem  there  probably  would  not  be  much  difficulty  identifying  the  first-mode  resonant 
frequency,  but  other  modes  might  be  hidden. 
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Reference  to  Shannon’s  sampling  theorem  shows  that  the  DFT  estimate  is  good  over 
the  frequency  band  of  0 to  l/(2At)  Hz  (Brigham,  1974).  It  is  important  to  note  that  the 
derivation  of  the  DFT  forces  both  the  time  and  frequency  domain  data  series  to  become 
periodic  with  a period  T,  even  if  this  is  not  the  case  in  reahty.  In  order  to  approximate  the 
infinite  length  of  signal  utilized  by  the  integral  transform,  the  discrete  transform  implicitly 
assumes  that  the  input  data  series  is  infinitely  repeated,  beginning  to  end.  If  the  beginning 
and  end  of  the  data  do  not  match  and  have  zero  slope,  an  infinitely  quick  jump  is  added  to 
the  data  series,  with  attendant  illegitimate  high  frequency  energy. 


2.3.2  The  DFT  and  Real-World  Systems 

As  presented  by  Aki  and  Richards  (1980),  a serious  problem  exists  when  estimating 
the  amplitude  spectrum  when  the  data  is  "contaminated"  by  noise.  The  DFT  model  of  the 
spectrum  does  not  expressly  take  the  presence  of  noise  into  account,  and  the  noise  must  be 
modeled  by  the  same  harmonic  sinusoids  as  the  signal  of  interest.  Since  real  field  data  is 
to  be  used,  noise  will  always  be  present  from  the  environment.  In  addition,  digitized  data 
will  always  have  "noise"  due  to  quantization  error.  The  Fourier  transform  of  the  noisy  data 
Nis 


N-FX 


(2.25) 


where  F = discrete  transform  of  the  signal 

X = discrete  transform  of  white  noise. 

Since  the  noise  is  a statistical  variable,  the  linear  combination  N is  also  a statistical 
variable.  Given  reasonable  approximations  (Bloomfield,  1976),  the  energy  spectral  density 
|Np  is  the  sum  of  two  squared  Gaussian  variables  (real  and  imaginary  parts  of  N)  and 
follows  the  chi-squared  distribution  with  two  degrees  of  freedom.  In  this  case  there  is  little 
confidence  that  the  sample  Fourier  spectrum  is  close  to  the  "true"  spectrum,  since  the 
standard  deviation  is  now  equal  to  the  mean.  There  is  little  certainty  whether  an  outlier  is 
a peak  value  of  site  amplification,  or  acceleration,  etc.,  or  a random  error. 

The  most  common  method  of  calculating  the  frequency  response  function  was  shown 
to  be  the  spectral  ratio,  Eq.  2.2.  A spectral  ratio  is  also  used  to  define  site  specific 
"amplification"  of  ground  motion  (Murphy  et  al.,  1971).  The  spectral  ratio  is  now  seen  to 
be  a ratio  of  two  chi-squared  variables  with  a common  mean  (from  a common  source)  and 
described  by  the  Fisher  F distribution  with  2x2  degrees  of  freedom  (Aki  and  Richards,  1980). 
There  is  very  little  statistical  certainty  with  so  few  degrees  of  freedom.  For  this  estimate 
there  is  a ninety  percent  probability  that  the  spectral  ratio  will  lie  between  0.053  and  19! 
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Simply  taking  more  data  is  of  no  help  in  reducing  uncertainty  since  the  only  effect  is 
to  increase  the  frequency  resolution  AF 


AF=^.  (2.26) 

NAt 

where  N = number  of  digitized  points 

At  = inverse  of  digitization  rate. 

Digitizing  at  a faster  rate  will  have  the  same  effect.  The  only  way  to  improve  the  certainty 
of  the  estimate  of  the  spectrum  is  to  average  adjacent  values,  so  that  the  variable  now  has 
four  degrees  of  freedom.  If  AF  is  very  fine,  many  values  can  be  averaged,  giving  a robust 
estimate  of  the  spectrum.  In  order  to  avoid  time-domain  aliasing,  the  number  of  frequency 
values  must  be  reduced  each  time  adjacent  frequency  bins  are  averaged.  The  trade-off  is 
confidence  for  frequency  resolution. 

When  an  average  of  adjacent  values  is  taken,  an  assumption  is  made  that  the  signal 
frequency  has  not  changed  between  the  two  bins.  Therefore,  if  the  signal  is  rapidly 
changing,  with  the  rate  of  change  being  on  the  order  of  the  desired  frequency  resolution, 
frequency  averaging  can  not  be  used  and  the  Fourier  spectrum  gives  a very  poor  estimate 
of  the  "actual"  spectrum.  Averaging  also  runs  opposite  to  the  need  to  maximize  frequency 
resolution.  This  conundrum  is  the  so-called  Uncertainty  Principle,  where  the  frequency 
resolution  is  inversely  proportional  to  the  length  of  time  signal,  which  in  turn  needs  to  be 
maximized  to  reduce  variance  (McClellan,  1982). 

A very  effective  approach  to  improving  the  statistical  reliability,  or  variance,  was 
introduced  by  Welch  and  employs  a form  of  ensemble  averaging  (Otnes  and  Enochson, 
1978;  Welch,  1967).  In  this  method  the  signal  is  broken  into  blocks,  the  spectral  estimate 
made  for  each  block,  and  the  resulting  spectra  averaged.  If  the  blocks  of  data  are 
overlapped,  e.g.  fifty  percent  when  the  Hamming  window  is  used,  the  bias  towards  the 
central  values  is  largely  counteracted.  This  method  can  only  be  implemented  if  there  is  a 
large  length  of  relatively  stationary  data  available. 

In  order  to  compare  and  understand  the  results  of  the  various  estimation  techniques 
discussed,  it  will  be  helpful  to  introduce  an  example  data  series  and  the  actual,  calculated 
spectrum  of  this  process.  The  time  history  is  shown  in  Fig.  2.7a,  and  is  taken  from  Kay  and 
Marple  (1981).  Figure  2.7b  shows  the  true,  calculated  spectrum  of  this  process.  The  process 
has  both  narrow-  and  wide-band  information,  with  two  line  components  very  close  together. 
While  the  spectrum  looks  simple,  this  combination  is  quite  difficult  to  estimate. 

Figure  2.8  shows  a typical  frequency  spectrum  calculated  using  the  periodogram 
method,  in  this  case  the  normalized  square  of  the  DFT.  The  oscillatory  nature  of  the 
estimate  is  a result  of  the  windowing  process.  Also,  leakage  causes  broadening  of  the  peaks 
which  can  merge  the  closely  adjacent  line  spectra,  shown  in  Fig.  2.7b. 
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Fig.  2.7a  The  example  data  process  to  test  spectral  estimates.  From  Kay  and  Marple  (1981). 
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Fig.  2.7b  The  true,  theoretical  spectrim  for  the  above  test  case  (Kay  and  Marple,  1981). 
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A much  improved  estimate  is  given  by  the  Welch  estimate  shown  in  Fig.  2.9.  In  this 
case  the  signal  has  been  windowed  so  that  the  large  side-lobes  (oscillatory  shapes  in  Fig.  2.8) 
have  been  reduced.  The  variance  has  been  greatly  reduced  by  the  averaging  process. 
However,  the  two  closely  spaced  line  spectra  expected  around  0.2  have  been  merged  due 
to  the  loss  of  resolution.  In  this  estimate  the  locations  of  the  spectral  peaks  are  correct  and 
the  shape  of  the  broadband  portion  is  similar  to  the  true  spectrum. 


2.3.3  An  Improved  Non-Parametric  Estimator 

An  extremely  powerful  non-parametric  spectral  estimator  for  almost  stationary  time 
series  has  been  proposed  by  Thomson  (1982).  This  is  the  so-called  multi-taper  method.  The 
method  is  derived  for  short  time  series  which  may  contain  line  spectra  as  weU  as  wide-band 
components.  The  method  is  rationally  derived,  as  opposed  to  the  ad  hoc  windowing  and 
filtering  of  the  classical  approach. 

In  the  classical  method  described  above,  the  data  is  windowed  to  try  to  control 
leakage  (bias),  Fourier  transformed,  and  smoothed  to  reduce  the  variance.  The  initial 
v^ndowing  increases  variance  and  weights  the  data  from  the  middle  of  the  time  series  much 
more  heavily  than  equally  valid  data  from  the  beginning  and  end  of  the  series.  The 
smoothing  (frequency  averaging)  is  only  rational  if  the  actual  spectrum  is  smooth.  Finally, 
the  information  from  the  phase  spectrum  is  discarded.  The  periodogram  estimates  are  not 
a "sufficient  statistic"  of  the  data  due  to  the  phase  information  being  discarded.  Any  finite 
spectral  estimate  is  an  under-determined  problem  (Thomson,  1982): 

Since  this  equation  {the  Fourier-based  spearum  - ed}  is  the  frequency-domain  expression 
of  the  projection  from  the  infinite  stationary  sequence  generated  by  the  random  orthogonal 
measure  dZ(f)  onto  the  finite  sample,  it  does  not  have  an  inverse;  hence  it  is  not  possible  to 
obtain  exact  or  unique  solutions.  What  we  desire  are  the  statistics  of  those  approximate 
solutions  that  are  both  statistically  and  numerically  plausible. 

Thomson  explicitly  deals  with  these  problems  by  proposing  a unified  method  which 
justifies  the  data  windows  used,  gives  consistent  estimates,  eliminates  bias  against  low 
amplitude  areas,  contains  separate  metrics  for  the  variance  of  line  and  broadband 
components,  and  for  which  the  size  of  the  time  series  enters  into  the  method  directly.  The 
solution  gives  local  independent  estimates  so  that  the  spectrum  at  a given  frequency  does 
not  depend  on  the  spectrum  at  a distant  frequency.  Both  the  spectrum  and  the  log  of  the 
spectrum  are  "good"  estimates,  while  for  the  periodogram  method  the  log  spectrum  is 
impossible  since  use  of  lag  windows  when  fi'equency  data  is  rapidly  changing  yields  negative 
estimate  values. 

The  multi-taper  method  starts  with  the  Cramer  spectral  representation  of  a function 
and  estimates  the  solution  to  this  integral  equation  by  a complex  orthonormal  eigenvector 
expansion.  The  user  chooses  the  number  of  expansion  terms  (tapers)  used:  the  more  terms 
used  the  less  biased  the  estimate,  at  the  expense  of  frequency  resolution.  The  class  of 
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Fig.  2.8  Periodogram  estimate  made  from  the  absolute  square  of  the  DFT.  The  typical 
oscillatory  shape  is  a result  of  windowing. 


Fig.  2.9  Classical  spectral  estimation  using  the  Welch  method.  The  spectrum 
is  to  be  compared  to  the  actual  spectrum  shown  in  Fig.  2.7b. 
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realizable  eigenvalues  is  limited  by  weighting  the  expansion  coefficients  by  prolate  spheroidal 
wave  functions,  which  happen  to  be  the  eigenfunctions  of  the  Dirichlet  kernel,  so  the  method 
stands  as  a coherent  whole. 

The  solution  is  reached  over  individual  local  intervals  surrounding  the 
eigenfrequencies  f,  i.e.  (f±w).  The  resulting  estimate  is  chi-square  distributed,  but  with  (2*# 
of  tapers)  degrees  of  freedom.  This  is  equivalent  to  (4Nw)  DOF,  where  N is  the  length  of 
time  series,  and  w is  the  frequency  radius  of  interest.  With  the  use  of  a reasonable  number 
of  tapers,  the  effective  window  is  very  rectangular  over  the  given  frequency  bin,  and  there 
are  very  low  side  lobes  (-80  dB  for  5 tapers)  effectively  eliminating  bias.  A simple  check  to 
see  if  too  many  tapers  are  being  used  is  to  examine  the  eigenvalues  in  full  precision;  if  they 
exceed  unity,  too  many  tapers  are  being  used. 

The  method  is  non-parametric  since  the  estimation  is  based  exclusively  on  the  time 
series  supplied  rather  than  on  a particular  model  of  the  process  producing  the  data.  If  any 
a priori  information  is  known  about  the  signal,  e.g.  whether  the  data  is  bandwidth-limited, 
or  that  no  line  spectra  are  present,  or  the  exact  nature  of  the  noise,  then  "estimations  of 
higher  apparent  resolution  can  be  made."  The  power  of  Thomson’s  method  is  that  very 
good  estimations  can  be  made  without  making  imprecise  a priori  estimations  that  can  skew 
the  results  towards  often  arbitrarily  predetermined  results.  One  of  the  major  problems  of 
traditional  analysis  is  that  the  researcher  processes  and  processes  the  data  until  it  yields  the 
results  that  were  expected.  However,  if  good  prior  information  about  the  system  is  available, 
a method  optimized  for  the  particular  situation  should  give  results  superior  to  those  from 
a general  solution. 

The  power  of  this  method  for  a knowledgeable  user  is  shown  in  Fig.  2.10,  which  is 
from  Thomson’s  paper  (1982).  The  estimate  for  the  test  time  series  is  for  all  practical 
purposes  the  same  as  the  true  spectrum.  However,  Thomson  did  much  more  than  just 
apply  his  algorithm,  which  gives  the  poor  result  shown  in  Fig.  2.11.  Thomson  utilized  the 
following  multi-step  estimation  procedure. 

(1)  Apply  the  multi- taper  algorithm. 

(2)  Calculate  the  variance  using  the  F-statistic  and  identify  the  frequencies  of  suspected 
line  components.  This  step  takes  a fair  amount  of  experience  and  knowledge. 
Thomson  makes  decisions  that  a user  of  lesser  experience  might  not  be  able  to  make. 

(3)  Subtract  the  effects  of  the  line  components. 

(4)  Further  pre-whiten  the  spectrum  with  an  autoregressive  (AR)  prediction  error  filter. 
Thomson  uses  a fifth-order  filter.  Pre-whitening  is  a method  used  to  remove  bias 
from  the  estimate  (Hardin,  1986;  Newland,  1984).  By  removing  spectral  peaks,  bias 
due  to  leakage  is  reduced  since  bias  is  proportional  to  the  second  derivative  of  the 
spectrum.  In  addition,  since  the  variance  of  the  estimate  is  proportional  to  the  mean 
of  the  estimate,  making  the  estimate  smaller  reduces  the  error.  By  making  the 
spectral  estimate  flat  — white  — bias  is  reduced  towards  zero  and  the  variance  is 
constant  for  all  frequencies.  An  initial  spectral  estimate  is  made  to  locate  peak 
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Fig.  2.10  Spectral  estimation  of  the  test  data  made  by  Thomson  using  his  multi-taper  method 
and  preprocessing  (Thomson,  1982). 


Fig.  2.1 1 Spectral  estimate  of  the  test  data  made  by  applying  multi-taper  method  with  no  preprocessing. 
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frequencies,  a filter  designed  to  remove  these  peaks  (prediction  error  filter),  new 
spectral  estimate  made,  and  the  estimate  postdarkened  by  recombination  of  the 
removed  peaks. 

(5)  Residual  spectrum  postdarkened  and  the  line  components  added,  giving  the  final 
estimate. 


2.4  Parametric  Methods  of  Process  Estimation 
2.4.1  Introduction 

For  the  non-parametric  approach  to  spectral  estimation,  no  assumption  was  made 
about  the  nature  of  the  data  except  that  the  time  series  values  are  identically  zero  outside 
the  windowed  section.  This  condition  is  true  for  a complete  transient,  but  obviously  is  not 
for  a stationary  segment  of  a transient.  Very  often  some  information  is  known  about  the 
signal  or  the  source.  It  might  be  something  as  simple  as  the  fact  that  the  frequency  content 
of  the  region  of  interest  is  band-limited.  For  example,  the  resonant  frequency  of  a building 
is  known  to  lie  within  a narrow  frequency  range.  Use  of  a priori  information  can  allow  a 
very  appropriate  model  of  the  process  to  be  used  to  estimate  a very  accurate  spectrum  with 
a small  amount  of  data. 

The  energy  spectral  density,  as  represented  by  the  standard  periodogram  approach 
discussed  above,  can  be  shown  to  be  identical  to  a parametric  model  of  a least  squares  fit 
of  the  time  series  to  a simple  harmonic  model  — the  DFT  (Kay  and  Marple,  1981).  The 
discrete  Fourier  frequencies  are  preassigned,  as  well  as  the  number  of  frequency  bins,  based 
on  the  digitized  period  and  digitization  speed  used  to  sample  the  data.  In  addition,  noise 
is  not  accounted  for  in  the  model,  the  energy  of  which  is  included  in  the  frequency  estimates. 
The  effects  of  noise  must  be  removed  through  the  various  averaging  schemes  discussed. 

In  the  above  discussion  of  Thomson’s  approach,  it  was  pointed  out  that  an 
unresolvable  limitation  of  the  classical  DFT  approach  to  spectral  estimation  is  the  fact  that 
a finite  set  of  values  (the  frequency  domain  time  series)  and  observations  are  used  to 
represent  a function  in  actuality  continuous  in  both  the  time  and  frequency  domains.  In  the 
parametric  approach,  a model  with  a finite  number  of  parameters  characterizes  the  process. 
The  recorded  data  is  used  to  estimate  the  parameters  of  the  chosen  model.  Note  that  there 
is  an  implicit  requirement  that  the  model  be  a good  representation  of  the  actual  physical 
process  being  studied. 
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2.4.2  Autoregressive  (AR)  Models 


The  most  common  approach  is  to  model  the  system  as  a linear  filter,  given  by  Eq. 

2.27 


V,-"  V.-i 


(2.27) 


where 


yj 

t 


= actual  data  sequence  (modeled  as  the  filter  output) 

= system  input  sequence  (assume  white  noise  for  spectral  estimation) 
= time  step  counter. 


The  output  is  seen  as  a combination  of  the  input  white  noise  history  acted  upon  by  the  "b" 
coefficients  plus  the  past  outputs  acted  upon  by  the  "a"  coefficients.  The  input  series, 
involving  the  "b"  coefficients,  is  a causal  moving  average  (MA)  process  (convolutional).  The 
series  involving  weighted  past  output  values  ("a"  coefficients)  is  a noncausal  autoregressive 
(AR)  process.  The  lengths  of  the  MA  and  AR  processes  must  be  explicitly  chosen  so  that 
the  model  best  represents  the  actual  process  (an  additional  piece  of  required  information 
compared  to  the  DFT). 

In  order  to  solve  for  the  model  parameters,  Eq.  2.27  is  rearranged  to  give 

9 p 

k=0  jfc=l 


The  Fourier  transform  of  Eq.  2.28  is  taken  by  applying  the  shifting  theorem  to  yield 

Applying  the  Z-transform,  where  2*  = e^,  to  Eq.  2.29,  and  rearranging,  gives 

„ (2.30) 


The  amplitude  spectrum  of  the  actual  data  now  becomes 


(2.31) 


since  the  amplitude  spectrum  of  the  white  noise  input  is  a constant  equal  to  one.  The 
numerator  polynomial  (in  order  q)  is  the  MA  process  while  the  denominator  polynomial 
contains  the  AR  coefficients. 
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If  a process  is  thought  to  be  a function  only  of  the  white  noise  input  being  acted  upon 
by  the  filter  model  the  process  is  said  to  be  a moving  average  process.  The  denominator 
of  Eq.  2.31  becomes  identically  one,  and  the  roots  of  the  factored  numerator  (MA) 
polynomial  are  all  "zeros"  of  the  equation.  If  the  roots  are  close  to  the  unit  circle,  then  the 
amplitudes  of  the  frequencies  close  to  the  zeros  will  become  very  small.  The  spectrum  will 
be  able  to  model  smooth  spectra  with  sharp  notches  very  well,  but  will  poorly  define  sharp 
peaks  unless  a very  large  number  of  parameters  are  solved  for.  The  AR  model  is  popular 
since  it  can  describe  a complex  process  with  very  few  parameters  calculated  from  a small 
length  of  data,  and  there  have  been  many  improvements  of  the  estimation  procedure. 

The  standard  periodogram  approach  to  spectral  estimation  can  be  shown  to  be  a 
special  case  of  the  MA  process  (Cadzow,  1982).  This  is  reflected  in  the  fact  that  for  the 
time  domain  analog,  the  autocorrelation  (from  the  Weiner-Khintchine  theorem),  it  is  only 
possible  to  calculate  a finite  number  of  time  lags  (q+1)  for  finite  data.  The  autocorrelation 
was  assumed  to  be  zero  outside  the  time  of  interest,  with  attendant  leakage  problems.  For 
broad-band  signals,  this  assumption  is  "almost  correct",  and  it  was  just  shown  that  a MA 
estimation  would  work  very  well  for  such  a process. 

Since  any  parametric  approach  is  useful  only  if  it  actually  is  a representative  model 
of  the  data  process,  it  is  important  to  develop  a feel  for  AR  processes.  It  has  been  shown 
above  that  an  AR  process  is  a feedback  process,  described  in  the  time  domain  as 


(2.32) 


The  current  output  of  the  model  is  a function  of  the  current  input  (assumed  to  be  random 
for  spectral  estimation)  and  a weighted  sum  of  the  past  outputs.  In  this  sense,  the  AR  filter 
can  be  seen  to  model  discrete  integration,  combining  past  outputs,  while  the  MA  filter 
models  a difference  (Robinson,  1982). 

Robinson  (1982)  describes  the  derivation  Yule  used  in  1927  (Yule,  1927)  to  formulate 
the  AR  spectrum.  Yule  imagined  a simple,  damped  SDOF  oscillator  — in  this  case  an  air- 
damped  pendulum.  This  system  can  be  described  in  a discrete  sense  by  the  homogeneous 
difference  equation 


c(f)  +a^c{t  - 1)  =0 


(2.33) 


where  c(t) 

Si 


displacement  amplitude  at  time  t 

model  parameters,  a2  is  the  reflection  coefficient. 
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(2.34) 


The  solution  to  Eq.  2.33  can  be  shown  to  be 


sin(f+l)  cOq 

c=e  

sin  (Oq 


where 


(2.35) 


CDQ=tan”^ 


(2.36) 


and  Cj  = impulse  response  function 

(j)o  = fundamental  frequency  of  the  impulse  response  function. 

This  result  should  be  compared  to  Eq.  2.3,  where  an  equivalent  system  was  modeled  from 
mechanical  principles  rather  than  by  a parametric  feedback  process. 

The  pendulum  is  then  excited  by  a white  noise  driving  function.  Yule  postulated 
several  boys  with  pea  shooters  irregularly  pelting  the  pendulum  from  different  locations. 
The  graph  of  the  motion  of  the  pendulum  through  time  will  be  quite  smooth,  with  the 
amplitude  and  phase  varying  continuously  as  given  by 

y(n)+a^y(n-l)+a^n-2)=x(n)  (2.37) 


where  Xj  = white  noise  input  excitation. 

The  solution  for  Eq.  2.37  is 


N 

y(t)=Ylc(k)x(t-k)  (2-38) 

ifc=0 

where  N = number  of  measurements  of  amplitude  taken 

c(t)  = impulse  response  function  given  by  Eq.  2.34. 

The  model  uses  a finite  number  of  chosen  parameters,  aj  and  a2.  Using  the  data 
series  y(t),  the  coefficients  are  found  by  regressing  yj  on  the  past  of  yj  — self-regression  or 
autoregression.  The  solution  utilizes  the  least  squares  normal  equations  filled  with  the 
empirical  autocorrelation  values  for  the  data  series,  and  are  called  the  Yule-Walker 
equations.  They  result  in  a Toeplitz  matrix  and  the  equations  are  rapidly  solved  using 
Levinson  recursion  (Levinson,  1947). 
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The  AR  coefficients  can  now  be  used  to  estimate  both  the  amplitude  and  phase 
spectra  of  the  data  series.  The  AR  coefficients  model  the  process  in  the  time  domain  - they 
model  the  time  series.  The  frequency  response  of  the  model,  i.e.  filter,  can  be  realized  by 
transforming  the  model  parameters.  Evaluating  the  time  domain  results  around  the  unit 
circle,  i.e.  taking  the  Z-transform,  yields 


N 


C^=j:C(n)z"  = 


n=0 


l+a^z+a^z^ 


(2.39) 


The  AR  model  is  much  different  than  the  MA-type  approach  made  by  the 
periodogram,  which  assumed  a sinusoidal  process  with  added  white  noise.  If  the  amplitude 
of  the  white  noise  increases,  the  graph  will  look  confusing,  but  the  periodogram  will  do  a 
fairly  good  job  (given  enough  observations)  of  picking  out  the  sinusoid.  However,  the 
improved  resolution  given  by  the  AR  approach  is  dependant  on  the  signal-to-noise  ratio 
between  the  sinusoid  and  background  white  noise,  doing  a poor  job  in  noisy  environments 
(Marple  and  Lawrence,  1987) 

For  the  AR  case,  the  roots  of  the  factored  AR  polynomial  are  all  "poles."  For 
frequencies  adjacent  to  poles  near  the  unit  circle,  amplitudes  will  be  very  large,  reminiscent 
of  the  shape  of  a circus  tent  close  to  the  poles.  If  it  is  expected  that  the  spectrum  is 
dominated  by  sharp  spikes,  then  the  process  can  be  well  modeled  as  an  autoregressive 
process. 

A serious  problem  can  occur  if  the  actual  autocorrelation  is  not  zero  outside  the 
limited  number  of  lags  available.  This  problem  occurs  when  the  process  is  actually  made 
up  of  sinusoid  and  white  noise.  A more  suitable  model  for  this  system  would  be  one  that 
does  not  window  the  signal,  i.e.  does  not  violently  truncate  the  autocorrelation.  The  most 
common  such  model  is  the  AR  model.  This  model  is  merely  Eq.  2.31  with  the  numerator 
terms  are  set  to  zero,  except  for  the  zero  time  lag  which  is  unity  (from  autocorrelation  of 
the  input  white  noise  series  xj. 

Insight  into  the  above  discussion  is  given  in  Fig.  2.12  (Marple  and  Lawrence,  1987). 
Figure  2.12a  shows  the  actual,  complete  autocorrelation  for  a single  sinusoid  on  the  left,  and 
the  true  power  spectral  density  for  a single  sinusoid  on  the  right.  The  assumptions  of  the 
periodogram  method,  and  the  realities  of  data-limited  discrete  processing,  yield  the 
truncated,  windowed,  autocorrelation  shown  on  the  left  of  Fig.  2.12b.  The  matching  biased, 
inconsistent,  periodogram  estimate  is  shown  on  the  right.  Finally,  the  left  side  of  Fig.  2.12c 
shows  the  effective  extrapolation  to  the  autocorrelation  made  by  the  very  simple  Yule- 
Walker  solution  for  the  AR  coefficients.  This  extrapolation  results  in  the  sharp,  unbiased 
estimate  shown  on  the  right  side  of  Fig.  2.12c. 

Unfortunately,  there  are  an  infinite  number  of  valid  extrapolations  of  the  truncated 
autocorrelation.  The  Yule- Walker  approach  (Marple  and  Lawrence,  1987)  to  estimating  the 
AR  model  parameters  is  based  strictly  on  the  truncated  autocorrelation  and  does  not  extract 
the  maximum  information  out  of  the  available  data. 
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Fig.  2.12  Burg's  extrapolation  of  the  autocorrelation  series  (ACS):  (a)  actual  ACS  and  the  true  spectrum 

(b)  truncated  and  incomplete  ACS  due  to  limited  discrete  data  and  the  classical  estimation 

(c)  Burg's  ACS  extrapolation  and  improved  AR  spectral  estimation  (Marple,  1987). 
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2.4.3  Maximum  Entropy  Method 

The  Burg  algorithm  (Burg,  1975)  is  the  most  popular  method  for  estimating  an  AR 
process.  This  is  the  so-called  Maximum  Entropy  estimation.  As  for  all  AR  methods,  the 
spectral  estimation  process  is  two-fold.  The  first  step  is  the  calculation  of  the  model 
parameters  from  the  recorded  data.  This  is  done  in  the  time  domain,  and  is  where  the 
effective  extension  of  the  autocorrelation  manifests  itself.  The  second  step  is  to  transform 
the  AR  coefficients  into  the  spectral  estimate.  However,  the  Burg  approach  estimates  the 
AR  parameters  without  making  any  explicit  estimates  of  the  autocorrelation. 

In  the  time  domain,  there  are  only  a limited  number  of  lags  of  the  autocorrelation 
known,  and  an  infinite  number  of  possible  extrapolations.  The  problem  is  to  choose  the 
'"best"  representative  of  this  infinite  class.  Burg  chose  the  spectrum  which  corresponds  to 
the  most  random  time  series  whose  autocorrelation  matches  that  of  the  actual  data  series. 
Burg  generates  a new  data  set,  based  on  the  actual  data,  which  is  the  most  random  possible 
in  an  entropy  sense.  The  reasoning  is  that  this  imposes  the  fewest  possible  constraints  on 
the  solution,  minimizing  the  bias.  The  extrapolation  supplies  additional  information  so  that 
the  estimate  has  very  high  resolution  compared  to  other  methods,  and  is  optimal  for  short 
data  sets  (Hardin,  1986). 

In  this  case,  the  concept  of  entropy  has  to  do  with  measuring  the  information  content 
of  a "message",  or  combination  of  symbols  (Shannon,  1948).  Maximizing  information  is  the 
same  as  maximizing  choice  when  selecting  a message;  and  the  more  freedom,  the  more 
uncertainty  as  to  the  outcome  (Radoski  et  al.  1975).  The  entropy  of  a distribution  is 
maximum  when  all  probabilities  are  equal.  Since  Burg’s  new  data  set  has  the  most 
randomness  (maximum  entropy),  the  associated  spectral  estimate  has  a high  resolution  since 
it  uses  all  the  possible  information  contained  in  the  estimated  autocorrelation,  when 
extrapolating  beyond  the  limited  number  of  lags  to  the  "true"  autocorrelation. 

Radoski  states  that  the  entropy  function  to  be  maximized  is 

E = f Ln[S((j>)]do>  (2-40) 


where  E 

= entropy  function 

S(o)) 

= spectral  density 

O)o 

= Nyquist  frequency,  7r/NAt  radians 

At 

= digitization  interval. 

The  problem  to  be  solved  is  to  calculate  a spectrum  that  maximizes  the  entropy  function  E, 
subject  to  the  constraints  of  the  actual,  lag-limited,  autocorrelation  estimate.  For  an  infinite 
data  series,  the  Burg  estimation  and  the  periodogram  would  be  identical.  Rather  than  use 
the  (incomplete)  autocorrelation  estimates  to  calculate  the  AR  parameters,  i.e.  the  Yule- 
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Walker  equations,  Burg  uses  a least  squares  approach  to  minimize  the  forward  and 
backwards  prediction  error  with  respect  to  the  last  filter  element  a„,  and  solving  for  the 
remainder  of  the  parameters  by  Levinson  recursion.  A very  important  advantage  of  Burg’s 
technique  is  that  it  is  optimized  for  use  of  a very  small  number  of  data  points  to  yield  a 
robust  system  estimate. 

Since  the  introduction  of  the  Burg  method,  many  other  estimation  schemes  have  been 
proposed  to  improve  on  limitations.  Swingler  (1979)  points  out  that  the  bidirectionality 
constraint  limits  the  application  of  Burg’s  method.  Not  all  deterministic  signals,  such  as 
exponential  decays,  can  be  successfully  modeled.  This  failing  should  not  be  too  surprising 
since  Burg’s  derivation  was  for  stationary  signals.  Swingler  also  states  that  the  insistence  on 
using  Levinson  recursion  causes  small  frequency  shifts  for  sinusoid  under  certain  phase 
conditions.  Swingler  presents  a non-recursive  least  squares  estimate  that  is  functionally 
identical  to  that  proposed  by  Ulrych  and  Clayton  (1976). 

Spectra  estimated  by  the  AR  models  have  a typical  swept-peak  shape.  The  peaks  are 
located  at  the  correct  places,  but  their  shape  is  not  a function  of  the  signal  itself.  In  fact, 
the  number  of  peaks  calculated  is  a function  of  the  AR  order  chosen  by  the  user.  There  will 
be  approximately  one  peak  for  every  two  parameters  chosen,  which  follows  from  Yule’s 
derivation.  The  typical  AR  spectrum  shown  in  Fig.  2.13a  exhibits  the  "peaky"  shape  and 
good  estimation  of  sharp  peaks  typical  of  the  AR  estimate.  In  contrast.  Fig.  2.13b  shows  the 
sharp  drops  well  modeled  by  the  MA  spectral  estimate.  As  expected,  the  combined  AR-MA 
estimation  is  a combination  of  the  two.  Note  that  the  shape  of  the  peaks  is  a function  of 
the  estimation  procedure  chosen  and  not  merely  a function  of  the  data  itself. 

These  points  are  illustrated  in  Fig.  2.14,  which  shows  the  least  squares  AR  estimation 
of  the  spectra  presented  in  Fig.  2.7  (Marple  and  Lawrence,  1987).  This  estimate  does  a far 
superior  job  modeling  the  actual  spectrum  than  the  classical  approach,  but  excessive  peaks 
are  evident. 


2.4.4  Autoregressive-Moving  Average  (ARMA)  Model 

It  has  been  mentioned  above  that  the  AR  estimate  is  not  suitable  for  cases  with  a low 
signal-to-noise  ratio  (SNR),  since  the  all-pole  model  is  not  valid  for  "sinusoid  in  white  noise." 
While  most  strong-motion  records  have  a high  SNR,  an  important  insight  can  be  gained  by 
examining  this  problem  (Kay  and  Marple,  1981).  The  white  noise  corrupted  AR  process  is 
defined  as  : 


n 


(2.40) 


where  x^ 
Wn 


= pure  AR  process 

= observed  white  noise  vdth  mean  = 0 and  variance  aj'. 
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(b) 


Fig.  2.13  Typical  parametric  spectra  from  the  (a)  AR(4)  model,  (2)  MA(4)  model, 
and  (3)  ARMA(4,4)  model  (Marple,  1987). 


Fig.  2.14  AR  spectral  estimate  of  test  series  to  be  compared  with  the  true  spectrum  in  Fig.  7b. 
Note  the  typical  "peaky"  shape. 
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Since  the  white  noise  is  uncorrelated  with  the  process  Xj,,  it  can  be  shown  that 


PyiZ)  = 


o'^+oj'  Aiz)A* 

At 

lill 

(2.42) 


where  Py(z) 

A(z) 

* 

At 

ji 


= spectral  density  of  the  corrupted  AR  process  in  the  Z-domain 
= filter  parameters 
= complex  conjugation 
= digitization  rate. 

= variance  of  the  output. 


The  noteworthy  property  of  Eq.  2.42  is  the  fact  that  the  spectral  density  is 
characterized  by  parameters  in  both  the  numerator  and  denominator.  This  model  was 
described  by  Eq.  2.31  and  is  an  ARMA  model  characterized  by  both  poles  and  zeroes.  The 
ARMA  model  is  very  powerful  in  that  it  can  easily  model  sharp  drops,  sharp  peaks,  and 
smooth  spectral  behavior.  It  is  also  the  most  parsimonious  estimator.  Since  the  ARMA 
model  is  the  most  general,  it’s  use  would  eliminate  the  problem  of  deciding  whether  the 
process  is  AR  or  MA. 


The  major  problem  of  implementing  an  ARMA  model  is  the  difficulty  in  calculating 
the  parameters.  The  most  straight-forward  approach  is  to  solve  the  Yule-Walker  normal 
equations  for  the  ARMA  process.  However,  in  this  case  the  equations  are  very  nonlinear 
and  difficult  to  solve  (Cadzow,  1982).  Simply  applying  a least  squares  solution  is  not 
sufficient  since  there  is  no  guarantee  of  convergence,  or  if  the  given  answer  is  the  result  of 
convergence  on  a local  extreme. 


A common  solution  is  to  solve  for  the  AR  and  MA  parameters  separately  and  then 
rationally  combine  the  two.  This  is  done  by  generating  the  extended  Yule- Walker  equations, 
an  example  of  which  is  given  in  Cadzow  (1982).  Treitel  derived  a least  squares  technique 
that  solves  the  ARMA  parameters  in  a unified  manner,  has  a minimum-delay  denominator, 
and  always  converges  (Gutowski  et  al.,  1978;  Treitel  et  al.,  1977). 

Another  problem  to  be  considered  is  the  estimation  of  orders  of  both  the  AR  and 
MA  parts  of  the  ARMA  model.  In  addition,  for  the  more  complete  system  identification 
models,  the  orders  of  the  noise  system  must  be  consistently  estimated  independently  from 
the  system  (van  den  Boom  and  van  den  Enden,  1974).  For  some  physical  systems,  such  as 
the  SDOF  oscillator  for  which  the  AR  model  was  derived,  the  estimation  of  AR  order  has 
physical  meaning.  For  other  processes,  the  order  must  be  guessed  at,  with  various 
theoretical  and  empirical  methods  of  deciding  the  optimum  AR  order.  Akaike  (e.g.  1981, 
1970)  approaches  the  problem  from  a similar  point  of  view  as  Burg,  and  uses  a maximum 
entropy  criteria  for  deciding  the  proper  model  order.  This  criteria  minimizes  the 
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"information  distance"  between  the  model  and  the  actual  process.  Too  few  parameters  and 
the  distance  is  great,  too  many  and  the  distance  slowly  increases  due  to  over-determination. 
The  criterion  is  the  so-called  Akaike  Information  Theoretic  Criterion,  or  AIC. 

There  are  many  empirical  methods  proposed  to  determine  the  proper  model  orders 
(e.g.  Marple  and  Lawrence,  1987;  Bohlin,  1984;  Astrom,  1980;  van  den  Boom  and  van  den 
Enden,  1974).  The  most  straight-forward  method  is  to  increase  the  orders  until  the 
"innovations"  series 


(2.43) 


where  yj 

5^ 


= actual  output  at  time  t 

= prediction  of  output  at  time  t made  at  time  t-1 


becomes  white  noise.  A second  check  that  all  the  information  is  being  "used"  by  the 
estimation  is  to  check  if  the  cross-correlation  of  the  input  series  and  the  innovations  is  white 
(Astrom,  1980).  At  this  point,  all  the  information  available  has  been  used. 

Through  a process  called  Wold  decomposition  it  can  be  shown  that  AR,  MA,  and 
ARMA  processes  are  related,  in  that  a ARMA  or  MA  process  of  limited  length  is  equivalent 
to  an  AR  model  of  infinite  length.  An  infinite  length  of  MA  filter  is  also  equivalent  to  a 
given  AR  or  ARMA  series.  Therefore,  a AR  model  with  order  much  greater  than  the  true 
model  order  is  a possible  solution  to  a low  SNR  condition.  In  order  to  avoid  the  problem 
of  spurious  peaks  in  the  spectral  estimate,  the  model  order  should  be  kept  to  less  than  half 
the  number  of  data  points. 

Another  method  for  dealing  with  noisy  data  is  to  compensate  for  the  uncorrelated 
noise.  This  is  done  by  either  weighting  the  zero  lag  of  the  autocorrelation,  or  a similar 
process  with  the  AR  coefficients.  The  rationale  behind  the  increased  weighting  of  the  zero 
lag  is  that  the  noise  is  considered  white,  with  an  autocorrelation  of  one  at  the  zero  lag  and 
zero  thereafter.  Increasing  the  value  of  the  actual  autocorrelation  at  lag  zero  effectively 
suppresses  the  effect  of  the  noise  on  the  solution. 

2.5  Spectral  Estimation 
2.5.1  Introduction 

The  parametric  and  non-parametric  methods  discussed  are  just  a few  of  the  more 
popular  models  for  describing  a process.  While  the  thrust  has  been  towards  spectral 
estimation,  it  must  be  remembered  that  all  these  methods  are  time  domain  solutions  with 
spectral  analogs.  The  important  question  now  becomes,  in  the  words  of  John  Tukey,  "When 
should  which  spectrum  approach  be  used?"  (Tukey,  1984).  Unfortunately,  there  is  no  easy 
answer  to  this  question  for  any  interesting  real-life  situations.  For  problems  of  interest,  such 
as  analysis  of  an  earthquake  strong  motion  record,  the  theoretical  caveats  needed  for  the 
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mathematical  derivation  are  not  applicable,  in  particular  the  assumption  of  stationarity. 

An  important  argument  is  made  by  Brillinger  and  Tukey  to  approach  the  subject  by 
using  "leading  cases"  (Brillinger  and  Tukey,  1984).  An  example  of  this  concept  is  the  use  of 
point  mass  assumptions  in  mechanics,  which  strictly  is  not  "correct,"  but  is  commonly  used 
without  complaint  and  with  outstanding  results.  Rather  than  approach  spectral  analysis  with 
strict  hypotheses,  Brillinger  and  Tukey  suggest  looking  at  the  practice  of  stochastic  spectral 
estimation  as  "...an  umbra  within  a penumbra." 

Three  successively  larger  spheres  of  application  can  be  described: 

(1)  An  inner  core  of  mathematically  derived  application.  This  region  is  generally  too 
narrow  for  application  to  real  problems.  An  example  is  a process  that  actually  is  a 
sum  of  a few  pure  sinusoids,  or  a process  with  these  qualities  unvarying,  realization 
after  realization.  No  example  from  physical  science  can  be  given. 

(2)  A middle  area  where  there  is  reasonable  understanding  of  the  process  and  the 
performance  of  the  method  used  to  assess  the  process;  the  meaning  of  the  spectrum 
is  clear.  There  usually  is  additional  information  about  the  process  at  hand,  but  all 
detailed  questions  are  not  yet  able  to  be  addressed.  This  large  area  is  dominated  by 
"stationary"  processes  with  finite  variance.  "This  is  not  because  the  region  of  use  of 
the  concept  of  a spectrum  is  confined  to  stationary  processes,  but  because  stationary 
processes  are  easier  to  think  about." 

(3)  The  third,  largest,  region  is  where  the  technique  is  actively  being  used  by  practitioners 
and  researchers.  In  this  area  the  understanding  of  the  model  and  process  gained  in 
the  inner  two  areas  is  applied  and  used.  A prime  example  is  the  use  of  spectral 
techniques  with  non-stationary  data  or  "unique",  one-of-a-kind  data  from  a process 
that  "...having  a process  would  require  an  ensemble  of  parallel  worlds." 
(earthquakes?).  The  famous  work  of  Munk  on  ocean  swells  is  a perfect  example  of 
the  important  information  that  can  be  gained  by  applying  spectral  methods  to  very 
non-stationary  time  series  (Munk  and  Snodgrass,  1957). 

A "non-stationary"  process,  such  as  an  earthquake  strong  motion  record,  can  be 
interpreted  in  many  different  ways.  On  one  end  of  the  scale  the  signal  can  be  seen  to  be 
made  up  of  short  stationary  segments,  each  of  which  has  a story  to  tell.  In  fact,  this  is  a 
common  and  useful  method  for  working  with  non-stationary  signals.  Both  the  Thomson  and 
Burg  estimates  of  spectral  density  were  derived  for  use  on  such  short  data  series.  Another 
approach  to  the  "microzonation"  of  a data  record  is  the  concept  of  evolutionary  spectra  put 
forth  by  Priestly  (1988).  In  this  approach,  the  spectrum  of  an  oscillatory  process  is  estimated 
in  the  "neighborhood  of  time  instant  t"  rather  than  over  all  time. 
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The  other  end  of  the  scale  is  to  view  physical  processes  in  a very  macroscopic  light 
(Brillinger  and  Tukey,  1984): 

All  we  know  of  the  world  is  consistent  with  the  idea  that  all  events  are  periodic  with  period 
10^®  years.  And  a process  made  up  of  all  displacements  of  a periodic  phenomenon  — with 
uniform  probability  — is  a stationary  process.  Thus  anything  we  find  in  the  world  could,  with 
this  definition,  come  from  a stationary  process. 


2.5.2  The  Concept  of  a "Spectrum" 

Non-stationarity  does  not  doom  spectral  analysis,  but  requires  extra  diligence  and  care 
to  insure  that  meaningful  estimates  are  given.  The  implications  of  each  step  must  be 
examined  so  that  effects  such  as  leakage  do  not  cover  the  desired  information  in  the  process. 
The  most  important  implication  that  must  be  addressed  is,  what  is  meant  by  the  concept  of 
a "spectrum"? 

As  commonly  used,  the  frequency  spectrum  represents  a given  signal  in  time  as 
sinusoids  of  various  frequencies.  The  various  frequency  components  present  in  the  signal 
are  presented  in  their  relative  energy  amplitudes.  This  usage  implies  an  oscillatory  signal 
made  up  of  "discrete"  sinusoidal  components.  An  example  might  be  the  several  modes  of 
excitement  of  a simple  oscillator,  or  the  characteristic  frequencies  of  a "good"  rotating 
bearing  verses  "bad".  The  calculation  of  the  spectrum  is  merely  a mathematical  transform 
which  sometimes  makes  a desired  facet  of  the  data  easier  to  find.  In  this  sense,  the 
spectrum  is  no  different  than  the  use  of  semi-log  or  square-root  plots. 

This  concept  of  frequency  spectrum  is  in  direct  contrast  to  the  representation  of  an 
arbitrary  waveform  shape  by  its  Fourier  frequencies.  A single  square-wave  pulse  has  a well- 
known  spectrum,  the  sine  function,  made  up  of  theoretically  infinite  frequency  components. 
In  the  present  context  this  square-wave  does  not  have  "frequency"  since  it  is  a single  pulse 
without  multiple  oscillations.  Multiple  occurrences  are  needed  before  one  can  speak  about 
the  frequency  of  an  occurrence  and  interpret  the  information  yielded  as  the  energy 
contained  in  the  various  frequencies  of  oscillation. 

For  practical  applications  even  this  straight-forward  "definition"  of  the  spectrum  has 
conflicting  implications.  Priestly  (1967)  gives  the  example  of  an  exponentially  damped  sine 
wave,  a very  common  physical  realization  : 

y^  = e sm(2Tc^t)  (2-44) 


where  f^  = 

An  example  of  this  signal  at  ^ = 20  Hz  and  a=8  is  given  in  Fig.  2.15a.  The  Fourier 
spectrum  pair  is  shown  in  Fig.  2.15b,  where  it  is  seen  that  the  spectrum,  while  having  a peak 
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Fig.  2.15a  Time  history  of  damped  sinusoid  f(t)=e  sin20(2Tct). 


-81 

Fig.  2.15b  Fourier  amplitude  spectrum  of  damped  sinusoid  f(t)=e  sin20(27rt). 
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at  20Hz,  contains  components  of  all  frequencies.  This  is  a function  of  the  damping,  which 
makes  the  signal  non-stationary  since  the  variance  is  continuously  changing.  However, 
practical  considerations  gives  the  alternative  interpretation  that  the  frequency  is  a constant 
20  Hz  while  the  amplitude  is  changing  as  a function  of  time.  No  engineer  with  real-world 
experience  would  dispute  the  pronouncement  that  the  signal  in  Fig.  2.15a  has  a spectrum 
indicating  a single  oscillatory  frequency  of  20  Hz.  In  this  case  the  results  are  "theoretically" 
incorrect,  but  for  the  engineer  estimating  the  process  parameters  the  "error"  at  0 Hz  has  the 
benefit  of  yielding  the  system  damping. 


2.5.3  Non-stationary  Signals 

It  is  interesting  to  note  that  if  an  "infinite"  sinusoid  at  20  hz  is  sampled  over  a finite 
time,  and  windowed  to  minimize  leakage,  results  are  very  similar  to  Fig.  2.15b.  In  essence 
sampling  turns  a signal  into  a transient  and  windowing  makes  the  transient  non-stationaiy. 
An  earthquake  strong  motion  record  becomes  slightly  more  "proper"  when  it  is  thought  of 
as  a self-windowed  signal  (of  infinite  length  if  desired)  with  some  white  background  noise. 
In  this  case,  the  effect  of  the  self-"windowing"  is  useful  since  it  can  provide  an  insight  into 
the  system  damping. 

Inherent  in  the  transformation  of  a time  signal  to  the  frequency  domain  is  the 
averaging  of  the  signal  components  over  the  sampling  period  T.  A piece  of  time  is  frozen 
over  this  period  and  the  assumption  made  that  all  time  before  and  after  is  the  same,  i.e. 
repeated  forever.  The  energies  present  at  each  component  frequency  are  integrated  over 
the  entire  time  period  T.  The  difficulty  with  non-stationary  signals  is  that  these  energies  are 
changing  during  this  period.  If  the  frequencies  present  are  changing  over  this  time  window, 
the  resulting  estimation,  regardless  of  method  used,  will  be  a smeared  average  as  if  all  the 
frequencies  with  energy  were  active  throughout  the  entire  period. 

For  weakly  non-stationary  processes,  the  effect  over  a small  time  period  is 
unimportant.  If  needed,  the  signal  can  be  cut  into  relatively  stationary  sections  and  spectra 
found  using  methods  specially  designed  for  short  data  segments,  i.e.  Thomson’s  or  Burg’s 
method.  This  approach  can  be  optimized  using  an  Akaike  AIC  criterion  (Akaike,  1974; 
Gersch  and  Brotherton,  1982). 

If  newer,  recursive/adaptive  methods,  such  as  Kalman  filtering,  are  used  to  yield  a 
time-varying  parameterization  of  the  process,  the  idea  of  "frequency"  must  again  be 
examined.  The  concept  of  an  "instantaneous"  frequency  in  the  traditional  context  is 
meaningless.  There  must  be  some  time  period  T over  which  to  sum  the  energy  expended 
(work)  at  each  frequency.  The  work  done  is  averaged  over  T,  hence  power  spectrum 
(power  = work  done/time  taken). 

This  problem  is  largely  one  of  semantics  arising  from  the  imprecise  manner  of 
everyday  language.  In  its  most  common,  and  traditionally  correct  usage,  the  spectrum 
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(amplitude,  energy,  or  power)  is  the  infinite  sum  of  sinusoids  with  arbitrary  amplitude  and 
phase,  necessary  to  represent  a given  waveform.  In  this  usage  the  waveform  can  be 
stationary,  or  a non-stationary  transient,  oscillatory  or  simple.  In  this  usage  a particular 
duration  of  signal  is  also  not  required,  e.g.  the  Dirac  delta  function  being  represented  by  a 
sum  of  sinusoid  of  all  frequencies.  As  discussed  earlier,  this  approach  is  strictly  correct  only 
for  the  infinite  analog  case,  and  application  to  the  discrete,  digital  domain  has  attendant 
problems. 

In  the  case  at  hand,  the  problem  is  to  describe  or  parameterize  a process.  In 
actuality,  all  the  approaches  discussed  take  place  in  the  time  domain  (the  classical  Fourier 
approach  is  actually  a computationally  simpler  version  of  the  autocorrelation  approach), 
yielding  an  impulse  response  function.  This  function  is  customarily  transformed  into  the 
frequency  domain,  now  called  the  frequency  response  function  or  transfer  function,  since  this 
representation  is  visually  less  complicated  and  the  desired  system  characteristics  are  more 
obvious.  In  particular,  the  system  fundamental  frequency  is  exceptionally  obvious  whatever 
method  is  used  to  estimate  the  process  spectrum. 

The  Fourier-based  methods  attempt  to  characterize  the  process  autocorrelation  by 
assuming  a uniformly  spaced  sinusoidal  basis  function  with  no  explicit  parameters  involved. 
This  method  was  seen  to  work  fine  if  the  process  can  accurately  be  modeled  thusly.  The 
parametric  methods  attempt  to  explicitly  model  the  process  through  equations  involving  a 
limited  number  of  parameters  that  must  be  experimentally  discovered.  These  methods  work 
well  if  the  process  can  be  accurately  modeled  by  the  chosen  method,  with  the  added 
difficulty  of  determining  the  model  order.  In  both  of  these  methods  the  process  is 
characterized  by  the  autocorrelation  or  model  parameters.  The  spectral  presentation  is  only 
an  alternative  method  of  presenting  the  data  in  a more  palatable  manner. 


2.5.4  The  "Instantaneous"  Spectrum 

Against  this  background  it  is  sensible  to  speak  of  the  "instantaneous  spectrum,"  which 
now  means  a frequency  domain  representation  of  the  behavior  of  the  system  at  a given 
moment.  At  any  given  instant  a filter  (system)  can  be  said  to  behave  as,  for  example,  an  AR 
process  of  order  two,  with  the  two  parameters  representing  a given  damping  and  resonant 
frequency,  free  to  change  as  the  source  motion  changes  or  the  structure  undergoes  damage. 
As  described  in  Fig.  2.1,  what  is  estimated  is  a filter  that  converts  the  known  input  into  the 
measured  output.  The  different  methods  discussed  are  merely  different  schemes  for 
estimating  that  filter,  which  can  be  presented  in  the  time  domain  or  the  frequency  domain. 
With  this  clearer  understanding  in  mind,  methods  allowing  time  varying  system  parameters 
will  be  discussed. 

An  early  method  of  directly  analyzing  time  varying  signals  was  the  use  of  complex 
demodulation  (Brillinger,  1988;  Bloomfield,  1976).  For  a signal  with  slowly  changing 


40 


amplitude  and  phase,  complex  demodulation  will  supply  an  estimate  of  these  values.  The 
method  can  be  thought  of  as  harmonic  analysis  in  a local  time  span  t±At.  The  frequency 
of  interest,  which  must  be  identified  beforehand,  is  isolated  by  narrow  band-pass  filtering. 
The  signal  is  modulated  by  a complex  exponential  and  then  locally  smoothed  to  yield  the 
local  amplitude  at  the  frequency  in  question.  This  method,  which  is  rather  cumbersome  and 
involves  ad  hoc  filtering  and  knowledge  of  spectral  peaks,  has  not  been  used  for 
identification  purposes  and  is  included  here  for  completeness. 

Using  traditional  methods  to  estimate  strongly  non-stationary  processes  can  yield  very 
unsatisfactory  results,  as  was  seen  above  in  the  damped  cosine  example.  Priestly  (1967) 
describes  early  attempts  to  define  an  "instantaneous"  spectral  estimate,  starting  in  1952 
(Page,  1952).  Priestly  puts  forth  an  evolutionary  (time  dependant)  power  spectrum  which 
he  describes  as  local  energy  distributions  over  frequency.  Priestly  computes  an  evolutionary 
spectrum  using  spectral  windows  with  various  parameters,  similar  to  the  periodogram 
approach,  and  therefor  suffers  the  same  difficulties.  Recursive  least  squares  methods  are 
now  available  which  estimate  AR  or  ARMA  parameters  describing  the  behavior  of  the 
process  at  each  time  step  (Marple  and  Lawrence,  1987). 


2.6  Modeling  Non-stationaiy  Processes  - Adaptive  Ffltering 
2.6.1  Introduction 

For  practical  stationary  problems  the  most  straightforward  method  is  in  essence  least 
squares  deconvolution,  or  calculation  of  the  Wiener  optimum  filter  (Kanasewitch,  1981; 
Wiener,  1964).  This  method  assumes  knowledge  of  the  autocorrelation  of  the  input  signal 
and  the  cross-correlation  of  the  input  and  output  signal.  The  system  parameters  are  solved 
for  directly  in  the  time  domain  and  the  defining  response  function  transformed  into  the 
frequency  response  function.  The  problems  involved  are  the  estimation  of  the  various 
covariance  functions  without  an  infinite  length  of  data,  and  these  problems  have  been 
discussed  in  previous  sections.  This  method  is  valid  if  the  signal  or  filter  does  not  change 
during  the  period  of  interest  and  could  be  applied  to  stationary-segmented  data. 

The  field  of  adaptive  filtering  was  formed  to  model  non-stationary  processes.  As  the 
statistics  of  the  signal  change  through  time,  the  filter  "adapts"  to  the  changing  variance  with 
new  parameters  that  reflect  the  structure  of  the  system  at  that  point.  The  predicted  value 
for  the  next  time  step  can  be  compared  with  the  actual  value,  and  the  difference 
(innovations,  Eq.  2.43)  will  give  a measure  of  how  well  the  filter  is  doing  its  job.  The  term 
"innovations"  is  used  because  this  information  is  new  information  that  can  not  be  predicted 
by  the  model  at  this  particular  step. 

Autoregressive  parameters  can  be  sequentially  estimated  so  that  the  parameters  are 
adaptive  to  the  changing  nature  of  the  process  (Marple  and  Lawrence,  1987).  The  AR 
parameters  are  updated  after  each  data  point,  tracking  slowly  non-stationary  signals.  A 
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forgetting  factor,  commonly  a damped  negative  exponential,  is  used  so  that  older  data  carries 
less  and  less  weight.  The  spectral  estimation  can  be  made  at  any  time  step  by  evaluating  the 
AR  parameters  around  the  unit  circle,  giving  the  spectral  description  of  the  behavior  of  the 
process  at  that  time.  The  most  effective  variant  of  the  recursive  least  squares  algorithm  was 
developed  by  Ljung  (Ljung  and  Soderstrom,  1983;  Ljung  et  al.,  1978;  Falconer  and  Ljung, 
1978). 

2.6.2  Kalman  Filters 

The  most  popular  direct  adaptive  filter,  or  process  model,  is  the  so-called  Kalman 
filter  (Kalman,  1960;  Kalman  and  Bucy,  1961).  Sorenson  (1970)  points  out  that  the  Kalman 
approach  is  a direct  descendant  of  Gauss’s  least  squares,  except  now  neither  the  signal  nor 
the  noise  model  must  be  stationary  — the  state  may  change  from  sample  point  to  sample 
point.  The  solution  is  recursive  and  is  presented  in  state-space  (Soderstrom  and  Stoica, 
1989)  which  uses  differential  equations  rather  than  integral  equations  to  represent  the 
model. 


Nau  and  Oliver  (1979)  state  that  the  Kalman  filter  is  based  on  a dynamic  AR  model 
defined  by  "two  concurrent  random  equations  of  motion": 


(2.45) 


the  AR(p)  equations  of  motion,  and  the  "motions"  of  the  parameters. 


(2.46) 


where  p 


at 

bt 


= number  of  prior  observations  utilized, 

= vector  of  p prior  data  observations  Xt.i,Xj.2,...pCt-p? 

= vector  of  p AR  parameters, 

= Gaussian  white  noise  with  0 mean  and  variance 
= Gaussian  white  noise  with  0 mean  and  covariance  matrix  Q. 


Equation  2.46  estimates  a value  of  0,  comprised  of  p previous  parameters,  through 
a random  walk  equation.  The  estimate  uses  the  weighted  p previous  data  points,  and  yields 
a new  observation  when  added  to  a new  noise  value. 

The  problem  is  to  use  Xj_j  to  filter  and  make  good  estimates  of  what  values  of  occurred, 
to  estimate  future  values  of  <f>^  from  Eq.  2.46  and  then  use  Eq.  2.45  to  forecast  future  Xj’s, 
repeating,  as  needed,  the  cycle  of  filtering  and  forecasting.  (Nau  and  Oliver,  1979). 

The  least  squares  solution  solves  the  equations  so  that  the  innovations  (Eq.  2.43)  — new, 
dynamic  information  that  cannot  be  predicted  — are  minimized  in  a least  squares  sense  each 
time  step 
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Kalman  noted  that  all  relevant  information  about  all  the  previous  data  is  contained 
in  the  posterior  covariance  matrix  of  the  parameter  distributions 


(2.47) 


where  = entire  history  of  the  process  up  to  and  including  time  t. 


The  estimate  is  an  explicit  function  of  only  the  data  history  and  the  previous  best  estimate 
of  the  state,  so  that  there  is  no  need  to  use  linear  regression  to  completely  re-estimate  the 
AR  parameters  after  each  new  time  step. 

Graupe  (1989)  discusses  the  extension  of  the  Kalman  filter  to  take  into  account 
colored  measurement  noise,  rather  than  strictly  Gaussian  white  noise.  This  is  possible  if  the 
noise  is  still  uncorrelated  with  the  signal.  This  Augmented  Kalman  filter  becomes  a ARMA 
representation,  with  the  colored  noise  represented  by  the  MA  component  in  a manner 
conceptually  similar  to  Eq.  2.42.  The  theory  behind  the  Kalman  filter  can  be  manipulated 
to  yield  the  system  parameters  for  the  case  where  there  is  no  a priori  information  about  the 
noise,  and  even  when  there  is  no  information  about  the  input  signal.  However,  these 
techniques  seem  designed  for  communications  problems  where  the  engineer  actually  has  at 
least  some  additional  conceptual  information  about  the  system  being  studied. 

The  so-called  extended  Kalman  filter  has  been  very  successfully  applied  to  non-linear 
estimation  problems  (Ljung,  1979;  Astrom  and  Eykhoff,  1971).  The  manner  of  application 
is  actually  straight-forward.  The  Kalman  model  is  constantly  updating  its  estimation  of  the 
dynamic  process  by  examining  the  innovations.  The  dynamics  can  be  due  to  a changing 
input  or  noise  process,  or  it  can  be  due  to  the  system  itself  changing.  The  effect  is  a 
linearization  between  single  time  steps,  but  if  the  system  is  changing  slowly  compared  to  the 
time  step  used,  the  linearization  is  "invisible"  and  the  non-linear  behavior  is  well  modeled. 
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CHAPTER  3 ESTIMATION  OF  SOIL  PARAMETERS  USING  SYSTEM 

IDENTIFICATION  TECHNIQUES 


3.1  A Framework  of  Understanding 

3.1.1  Introduction 

This  chapter  will  present  some  interesting  and  fruitful  approaches  to  applying  the 
system  identification  (SI)  methods  presented  in  the  previous  chapter.  Given  an  excitation 
input  and  the  system  output,  the  goal  is  to  identify  various  physical  parameters  that  control 
the  dynamic  behavior  of  the  soil  mass.  Typical  parameters  to  be  identified  are  soil  layer 
resonant  frequency  and  damping  ratio.  These  are  the  obvious  choices  since  it  has  been 
shown  that  the  ARMA  evaluation  techniques  were  derived  to  model  the  dynamic  process 
as  SDOF  oscillators.  Other  parameters  that  has  been  studied  are  the  "amplification"  of 
ground  motion  as  the  seismic  energy  moves  from  stiff  bedrock  to  the  much  more  compliant 
soil,  and  identification  and  simulation  of  strong  ground  motions. 

This  report  will  be  concerned  with  the  application  of  SI  to  the  problem  of  identifying 
soil  parameters.  It  will  not  attempt  to  cover  the  large  field  of  work  done  in  the  application 
of  SI  to  the  soil-structure  interaction  problem  (e.g.  many  papers  in  the  journal  Earthquake 
Engineering  and  Structural  Dynamics;  Ghanem  et  al.,  1991).  This  report  will  present  some 
interesting  examples  of  the  use  of  dynamic  models  to  characterize  in  situ  soil  behavior. 

The  work  to  date  in  seismic  SI  can  be  divided  into  two  large  groups  — continuous 
excitation  and  transient  excitation.  Examples  of  continuous  excitation  are  shakers  and 
ambient  excitation.  Earthquakes  and  seismic  testing  methods  are  examples  of  transient 
excitation.  In  relation  to  analysis  methods  available,  these  two  group  can  each  be  divided 
into  stationary  and  non-stationary  sources.  Excitation  of  the  soil  through  a servo-controlled 
shaker  is  certainly  a known,  stationary  force,  so  that  the  straight-forward  methods  of  Bendat 
and  Piersol  (e.g.  1986)  can  be  used  to  calculate  soil  damping  and  resonant  frequency.  The 
use  of  ambient  vibrations,  while  often  considered  to  be  white  processes,  can  very  well  be  at 
least  weakly  non-stationary. 

Seismic  methods,  such  as  impulse  loading  and  use  of  explosives,  are  non-stationary 
but  very  repeatable.  The  repeatability  of  the  input  allows  ensemble  averaging,  and  the 
averaged  process  can  be  well  characterized  by  simple  Fourier  analysis.  Finally,  earthquake 
excitation  is  obviously  non-stationary  and  non-repeatable.  It  is  for  this  situation,  where  the 
amount  of  information  is  limited,  that  system  identification  techniques  are  of  the  most  use. 

When  analyzing  the  in  situ  soil  profile,  the  geotechnical  engineer  must  work  with  very 
large  "structures."  It  is  difficult  to  artificially  excite  a significant  mass  of  the  earth  to  levels 
suitable  to  employ  system  identification  techniques  for  strain  levels  even  approaching 
threshold  strain.  An  obvious  source  for  larger  strain  levels  is  an  earthquake.  In  the  past 
several  years  there  have  been  a few  instrumented  sites  where  strong-motion  records  have 


been  recorded  within  the  soil  as  well  as  on  the  surface.  For  these  cases  the  techniques 
discussed  in  the  previous  chapter  can  be  used  to  calculate  soil  parameters,  as  well  as  to  study 
how  the  parameters  might  change  during  the  shaking. 


3.2  Continuous  Forced  Excitation  of  the  Soil 
3.2.1  Introduction 

The  simplest  method  of  determining  the  dynamic  behavior  of  a structure  is  to  use 
forced  vibration.  In  this  case  the  excitation  is  under  complete  control  of  the  experimenter, 
and  the  input  forcing  function  known.  If  the  soil  can  be  modeled  as  a one  layer  system  for 
the  problem  at  hand,  then  solution  presented  in  Sect.  3.2.2  can  be  utilized  (Richart  et  al., 
1970).  The  driving  function  is  stationary,  and  duplicate  data  can  easily  be  collected.  This 
makes  a proper  DFT-based  spectral  estimate  very  accurate.  The  gain  factor  can  be  plotted, 
the  resonant  frequency  and  soil  stiffness  found,  and  the  soil  damping  calculated  from  the 
half-power  point  (Bendat  and  Piersol,  1980).  This  method  can  be  expanded  to  examine 
higher  modes  of  vibration. 


3.2.2  Applications  to  Multiple  Layered  Profiles 

In  the  more  common  situation  the  soil  must  be  considered  a layered  system,  with 
each  layer  having  its  own  distinct  set  of  descriptive  parameters.  This  case  is  modeled  as  a 
series  of  SDOF  linear  oscillators.  The  system  will  be  characterized  as  the  compound  filter 
that  converts  the  known  force  input  into  the  measured  surface  displacement  output.  The 
solution  to  this  inverse  problem  is  not  necessarily  unique,  especially  if  the  actual  input 
motion  is  not  precisely  known  and  near-by  "bedrock"  motion  is  used  (Udwadia  et  al.,  1978). 

An  N-layered  soil  system  typically  can  be  modeled  as  an  N-degree-of-freedom 
lumped  mass  system,  as  illustrated  in  Fig.  3.1.  If  the  mass  distribution  is  known,  Udwadia 
has  derived  a method  whereby  the  damping  and  stiffness  of  all  layers  can  be  uniquely 
determined  for  forced  vibration  with  one  surface  transducer  (1986).  The  method  is  valid  for 
linear  response  only.  Since  the  force-time  history  of  the  exciting  shaker  and  the  system 
output  is  known,  this  method  is  fundamentally  similar  to  the  simple  case  just  described. 
Udwadia  notes  that  a good  signal-to-noise  ratio  is  an  important  parameter  to  insure  low 
variance  in  the  system  parameter  estimates. 

A similar  but  less-controlled  variant  of  Udwadia’s  method  is  to  use  background 
ambient  vibrations  as  the  excitation,  making  the  potentially  incorrect  assumption  that  this 
excitation  is  gaussian  white  noise.  An  apparent  application  of  this  approach  is  presented  by 
Ohmachi  for  a suite  of  locations  in  the  San  Francisco  Bay  area  (Ohmachi  et  al.,  1991). 
Three-component  surface  measurements  were  taken,  three  times  each  at  the  several 
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Fig.  3.1  n-Degree-of-Freedom  lumped  mass  soil  system  (from  Udwadia,  1978). 


locations.  Small  subsections  (1024  points)  of  the  data  were  taken,  Fourier  transformed,  and 
smoothed  by  multiple  passes  of  a Hanning  lag  window  (see  Newland,  1984).  Three  such 
spectra  were  then  averaged  for  each  component  at  each  location. 

This  study  was  based  on  a method  developed  by  the  Japanese  National  Rail  (JNR) 
to  study  soil  characteristics  along  rail  lines  using  microtremors  (frequency  content  0.5  to  20 
Hz)  as  excitation  (Nakamura,  1989).  Noting  the  virtual  impossibility  of  taking  an  adequate 
amount  of  data  without  colored  interference  from  social  activities,  Nakamura  derives  a 
method  to  rid  the  data  of  these  effects.  He  claims  that  the  social  sources,  the  "noise", 
produce  mostly  Rayleigh  waves,  dominated  by  vertical  motion.  The  attempt  is  then  to 
eliminate  the  effects  of  the  surface  wave. 

Based  on  a comparison  of  five  rock  sites  and  four  soil  sites,  Nakamura  claims  that 
for  rock  the  magnitude  of  vertical  acceleration  is  the  same  as  the  horizontal.  For  soils  the 
horizontal  amplitude  is  greater  than  for  the  vertical.  This  is  a broad  assumption  based  on 
the  small  amount  of  data  presented.  He  goes  on  to  claim  that  the  amplification  of  motion 
from  a stiff  layer  to  soft  is  through  multiple  reflections  throughout  the  soil  layer.  This 
assumption  is  patently  false  since  it  violates  conservation  of  energy.  In  fact,  net  energy  is 
lost  through  reflection,  otherwise  an  amplifier  using  no  input  energy  could  be  built  with  a 
pile  of  plates. 

This  background  is  used  by  Nakamura  (1989)  to  state  that  the  effect  of  the  R-wave 
can  be  measured  by  taking  the  ratio  of  the  vertical  displacement  between  the  surface  and 
a subsurface  acceleration  record.  Breaking  the  previously  made  assumptions,  Nakamura 
proceeds  to  derive  a method  that  he  claims  will  allow  the  estimation  of  the  transfer  function 
of  "surface  layers  from  tremor  on  the  surface  only."  This  is  done  through  a ratio  of  the 
horizontal  and  vertical  spectra,  and  called  the  "spectral  amplification  factor."  The  trick 
seems  to  be  that,  for  the  frequencies  of  interest,  the  R-wave  has  the  same  spectral  amplitude 
in  both  the  vertical  and  horizontal  directions.  He  then  adds,  "Note,  however,  that  the 
estimation  accuracy  drops  when  there  exists  a noise  tremor  agreeing  with  the  prevailing 
frequency  in  the  estimated  transfer  function." 

While  this  technique  would  nullify  common  excitation  factors,  it  would  also  nullify  any 
common  response  characteristics.  Just  because  a white  excitation  itself  is  uncorrelated  (non- 
orthogonal)  does  not  mean  that  the  vertical  and  horizontal  system  response  will  be 
uncorrelated.  Until  further  explanation  is  published,  this  method  seems  to  be  biased  and 
not  connected  to  system  behavior  in  an  unique  manner.  Finn  (1991)  states,  "Nakamura’s 
procedure  for  determining  site  periods  and  amplification  factors  is  based  on  some  tenuous 
assumptions."  In  addition,  the  parameters  reported  — "predominant  frequencies"  and 
"amplification  factors"  — are  not  the  commonly  understood  process  parameters  and  are  index 
properties. 
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The  results  presented  in  Nakamura  (1989)  shows  that  the  method  gives  acceptable 
results  for  certain  conditions.  Indeed,  JNR  would  not  monitor  the  entire  length  of  all 
Shinkansen  lines  if  the  method  gave  no  useful  results.  Since  the  given  derivation  is 
contradictory,  either  the  method  works  because  of  R-wave  coincidentally  canceling  for  the 
sites  of  interest,  or  there  is  a more  complete  derivation  that  has  not  been  presented  in 
English.  In  addition,  there  is  a problem  with  using  microtremors  to  characterize  the  strong 
motion  behavior  of  a site.  Data  from  Japan  show  that  the  microtremor  techniques 
overestimate  the  transfer  function  by  up  to  a factor  of  five  compared  to  the  strong  motion 
results  (Seo,  1989).  In  addition,  evidence  shows  that  the  fundamental  period  shifts  for  strong 
motion  accelerations. 


3.2.3  Modal  Analysis 

The  behavior  of  a well-defined  structure  under  forced  vibration  can  be  used  to 
ascertain  information  about  the  supporting  soil  system.  The  results  of  these  tests  have  been 
used  to  measure  soil  stiffness  beneath  bridge  piers  (Crouse  and  Hushmand,  1987).  Thirty 
accelerometers  were  placed  at  various  locations  of  the  bridge,  and  a detailed  modal  analysis 
of  the  bridge  response  to  a large  eccentric  mass  shaker  was  undertaken.  It  was  found  that 
the  second  (primary  torsional)  mode  was  found  to  best  characterize  the  soil-structure 
interaction. 

A variational  form  of  Rayleigh’s  principle  was  used  to  estimate  the  stiffness  of  the 
Winkler  springs  used  to  model  the  stiffness  of  the  soil  beneath  the  bridge  footings. 
Lagrange’s  equations  were  used  to  write  the  equations  of  motion  relating  to  an  assumed 
deflection  shape.  The  actual  geometric,  inertial  and  displacement  responses  of  the  bridge 
were  then  substituted  into  these  equations,  and  spring  stiffnesses  estimated  when  the 
substitution  led  to  convergence. 

The  calculated  stiffness  of  the  sounding  soil  mass  matched  very  well  with  results  from 
finite  element  analysis  and  experimental  estimates.  The  computed  resonant  frequencies 
were  almost  identical  to  those  observed  for  the  first  four  modes.  While  yielding  good  results, 
this  method  is  quite  complicated,  both  experimentally  and  computationally.  In  this  case 
thirty  channels  of  acceleration  data  had  to  be  recorded  and  analyzed.  In  addition,  the 
authors  note  that  the  technique  will  only  work  with  very  well  characterized  structures  since 
the  effects  of  soil  and  structure  must  be  separated.  The  values  given  are  for  an  averaged 
lumped  soil  mass  and  not  for  each  soil  layer  independently,  therefore  this  method  would 
probably  not  be  of  much  help  to  identify  layers  of  potentially  liquefiable  sand. 


3.2.4  Estimation  of  Soil  Properties  from  Impedance  Functions 

One  powerful  method  of  describing  the  dynamic  interaction  between  a vibrating  rigid 
body  and  the  soil  is  through  dynamic  impedance  functions.  The  method  was  developed  to 
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analyze  rigid  machine  foundations  as  a function  of  the  excitation  frequency.  The  dynamic 
impedance  function,  K^,  is  decided  for  an  "associated"  massless  foundation,  and  then  the 
steady-state  behavior  of  the  actual  foundation  and  any  mass  resting  on  it  can  be  easily 
calculated.  The  impedance  function  can  be  determined  analytically,  numerically,  or 
experimentally.  The  following  discussion  of  impedance  functions  is  from  a paper  by  Gazetas 
(1983). 


For  each  harmonic  excitation  frequency  the  impedance,  K,  is  defined  as  the  ratio  of 
the  steady-state  force  and  resulting  displacement  at  the  interface  of  the  soil  and  foundation. 
The  vertical  impedance  can  be  written  as 


i(i>r 


Ve 


tat 


harmonic  vertical  force 
uniform  harmonic  settlement 


(3.1) 


Since  the  dynamic  force  and  displacement  will  generally  be  out  of  phase,  the  displacement 
can  be  split  into  one  component  which  is  in  phase  with  the  force,  and  an  orthogonal 
component.  The  dynamic  impedance  can  then  be  expressed  in  a complex  form: 

= (3.2) 

where  m = displacement  mode,  i.e.  vertical,  horizontal,  rotational,  torsional,  or  coupled. 

In  this  format,  the  real  part  contains  the  effects  of  the  stiffness  and  inertia  of  the  supporting 
soil  mass,  and  the  imaginary  part  describes  the  radiational  and  material  damping  of  the  soil 
mass.  Both  the  stiffness  and  material  damping  are  believed  to  be  frequency  independent, 
so  Ky  will  vary  due  to  the  frequency  dependant  behavior  of  inertial  mass  and  radiational 
energy  loss. 

The  SDOF  oscillator  can  also  be  represented  by  the  complex  dynamic  impedance 
function  in  the  form  of  Eq.  3.2: 

K((ji)=(k-M(ji^)  + iC(ji 

where  (k-Mo)^)  = 

Cg)  = K2 

M = mass 

k = static  soil  stiffness 

C = soil  damping. 

Therefore,  if  the  dynamic  impedance  function  for  a specific  footing  and  soil  system  can  be 
both  experimentally  measured  and  numerically  calculated,  the  soil  parameters  can  be 
identified  through  an  iterative  process. 
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The  impedance  function,  Eq.  3.3,  can  be  written  in  a more  common  form: 


(3.4) 


where  k 
13 


= static  soil  stiffriess 
= = C(Ojj/2k  = damping  ratio 

= y(k/M)  = natural  frequency 


G) 


n 


or,  in  short-hand  form. 


A'=ife:|K+z(i)c^| 


(3.5) 


where  k = static  soil  stiffness 

K = (l-o^/Wj,^)  = stiffness  coefficients 

Cj  = C/k.  = damping  coefficients. 


Gazetas  notes  that,  in  the  form  of  Eq.  3.5,  the  dynamic  impedance  of  the  SDOF  oscillator 
can  be  seen  as  the  product  of  the  static  stifbiess  of  the  soil  and  a complex  number 
incorporating  the  dynamic  characteristics  of  the  system.  The  stiffness  coefficients  decrease 
with  increasing  frequency  while  the  damping  coefficients  stay  constant,  as  shown  by  the  chart 
in  Fig.  3.2. 

In  reality,  however,  the  relation  between  the  stiffness  and  damping  coefficients,  and 
frequency,  is  complex  for  sod-foundation  systems.  The  behavior  of  the  dynamic  factor  of  the 
impedance  function  will  depend  on  the  mode  of  vibration,  geometry,  rigidity,  and 
embedment  of  the  foundation,  and  the  profile  and  elastic  properties  of  the  soil  system. 
These  complications  are  the  reason  that  difficult  numerical  techniques,  only  available  since 
the  mid-1970’s,  must  be  used  to  estimate  impedance  functions.  The  experimental 
determination  of  these  functions  has  also  been  a difficult  task  (Crouse  et  al.,  1990). 

A small  rigid  footing  can  be  excited  by  forced  vibrations  to  determine  actual 
impedance  functions  in  situ  (Luco  and  Wong,  1990).  The  soil  is  modeled  as  several 
horizontal  layers  overlying  a half-space.  The  parameters  of  the  theoretical  function  are 
varied  until  there  is  a least  squares  match  between  it  and  the  experimental  function.  This 
allows  the  identification  of  the  shear-wave  velocities  (stiffness),  material  damping  ratios,  and 
Poisson’s  ratios  of  the  soil  layers.  Luco  and  Wong  used  the  full  set  of  impedance  functions, 
incorporating  all  five  modes  of  vibration,  in  their  theoretical  calculations. 
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Fig.  3.2  Dynamic  stiffness  and  damping  coefficients  of  a SDOF  oscillator  (from  Gazetas, 
1983) 
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This  technique  was  evaluated  with  a set  of  data  from  a foundation  test  in  the  Imperial 
Valley,  CA  (Crouse  et  ah,  1990),  which  only  contained  the  experimental  functions  for  three 
modes.  The  resultant  velocities  were  different  from  the  experimentally  measured  values,  but 
within  the  experimental  variance  of  the  SASW  method  used  to  measure  the  in  situ  velocities. 
However,  the  damping  ratio  values  were  very  erratic  and  unrealistic. 

The  impedance  function  method  has  been  further  refined  to  improve  the  damping 
estimates  (Luco  and  Wong,  1992).  The  primary  concern  of  this  paper  was  to  ascertain  the 
effect  of  frequency  on  the  identification  accuracy.  Since  the  applied  force  is  proportional 
to  the  square  of  the  frequency,  the  force  at  low  frequencies  is  small  and  approaches  the 
noise  level.  However,  damping  values  are  best  defined  at  low  frequencies.  It  was  shown 
that  damping  ratios  for  soil  deeper  than  one-quarter  of  the  wavelength  associated  with  the 
frequency  of  interest  cannot  be  accurately  identified.  A similar  relationship  between 
wavelength/frequency  and  depth  of  sampling  was  found  in  the  derivation  of  the  SASW 
method. 

Luco  and  Wong  (1992)  also  found  that  data  error  degrades  the  estimates  of  the 
properties  of  deeper  layers,  especially  the  damping  ratio.  As  expected,  the  use  of  more 
impedance  functions  (different  modes)  improves  the  identification.  The  authors  improved 
their  results  by  the  use  of  a weighting  function  (smoothing)  on  successive  estimate  iterations. 

While  the  use  of  impedance  functions  to  estimate  in  situ  soil  properties  shows 
promise  and  has  good  theoretical  grounding,  it  is  very  limited  as  to  the  depth  of  soil  it  can 
characterize.  The  maximum  is  only  a few  meters.  The  experimental  and  analytic  overhead 
is  also  quite  daunting.  The  forced  vibration  techniques  differ  from  the  typical  system 
identification  methods  discussed  in  Chapter  2 in  that  the  parameters  were  identified  through 
mechanistic/deterministic  models  of  specific  behavior  rather  than  fitting  the  output  time 
series  to  a more  general  model  such  as  the  AR  feedback  model. 


3.3  Non-stationary  Excitation  of  the  Soil  System 
3.3.1  Continuous  Excitation 

Forced  vibration  methods  are  often  limited  because  there  is  no  way  to  input  enough 
energy  into  the  system  to  allow  a large  amount  of  soil  to  become  involved.  A solution  to 
this  problem  is  to  use  natural  excitation.  A constant,  repeatable  source  would  be 
appropriate,  since  it  would  allow  large  data  sets  and  redundancy  to  minimize  bias  and 
variance  in  the  estimates.  Ideally,  simple  DFT  methods  could  be  used  with  excellent  results. 
Of  course,  this  assumes  that  the  natural  signal  is  stationaiy,  as  was  assumed  by  Ohmachi  et 
al.  (1991).  For  example,  Gersch  et  al.  (1973)  used  a maximum  likelihood  ARMA  estimation 
of  the  system  parameters  from  two  different  length  time  series  from  a wind-loaded  building. 
The  estimated  damping  ratio  and  resonant  frequency,  and  their  coefficients  of  variation, 
were  consistent  for  the  first  two  modes  of  the  625  point  and  2500  point  series.  This  led  the 
authors  to  conclude  that,  in  this  case,  the  wind  loading  is  a white  noise  excitation. 


53 


Ambient  vibration,  however,  is  non-stationary  in  the  long-term.  A straight-forward 
solution  is  to  break  up  the  non-stationary  record  into  stationary  pieces.  The  simplest 
method  is  to  assume  that  the  variance  is  changing  slowly,  and  break  the  record  into  relatively 
short,  equal-length  segments.  These  segments  can  then  be  independently  analyzed  and 
modeled  (Popescu  and  Demetriu,  1990;  Gersch  and  Martinelli,  1979). 

A much  more  efficient  approach  is  to  use  the  Akaike  Information-theoretic  Criterion 
(AIC)  to  divide  the  record  into  locally  stationary,  contiguous  but  not  necessarily  equal  length 
segments  (Gersch  and  Brotherton,  1982).  Each  locally  stationary  segment  is  then  assumed 
to  be  made  up  of  two  components.  The  first  is  a stationary  time  series  that  is  the  same  for 
all  the  segments  and  represents  the  time  invariant  parameters  of  the  structure  being 
investigated.  The  second  is  a non-specific  auxiliary  signal  and  varies  from  segment  to 
segment,  and  can  be  due  to  various  noise  sources. 

While  it  was  shown  that  a regularly  sampled  vibration  record  is  exactly  characterized 
by  a 2n-2n  order  ARMA  model  (Gersch  et  al.,  1973),  Gersch  simplifies  the  problem  by 
noting  that  the  MA  parameters  contribute  very  little  to  the  estimates  of  the  structural 
parameters.  As  explained  in  Chapter  2,  the  problem  is  greatly  simplified  by  modeling  the 
process  as  strictly  AR,  since  the  calculation  procedure  is  well  characterized.  It  was  stated 
earlier  that  MA  parameters  provide  information  on  the  phase  relation  of  the  time  series, 
while  the  values  of  the  natural  frequency  and  damping  ratio  depend  only  on  the  AR 
parameters  (Chang  et  al.,  1982). 

The  model  used,  with  an  auxiliary  process,  is  the  ARX  model.  The  AR  spectral 
estimate  is  calculated  from  the  AR  parameters,  and  a corresponding  pole  plot  made  for 
each  optimal  time  segment.  Since  the  structural  part  is  common  to  all  segments  and  the 
auxiliary  (noise)  part  different  in  each,  the  complex  poles  of  the  structural  part  will  be  in  the 
same  location  while  those  of  the  noise  process  will  vary  in  each  pole  plot.  This  is  how  the 
parameters  of  interest  are  pulled  out  of  the  general  model. 

This  ingenious  method  allows  more  certain  estimates  of  the  damping  parameter. 
Since  the  variance  of  the  damping  estimate  is  a function  of  the  inverse  of  the  length  of  data 
used  to  make  the  estimate,  long  stretches  of  stationary  data  are  usually  required  for  robust 
estimates.  The  present  method  allows  the  estimates  from  each  segment  to  be  combined  with 
a final  variance  inversely  proportional  to  the  length  of  the  sum  of  segments. 

Gersch  and  Brotherton  (1982)  present  an  example  based  on  ambient  vibration  taken 
from  the  top  story  of  the  Jet  Propulsion  Laboratory  at  the  California  Institute  of  Technology. 
The  values  for  natural  frequency  and  damping  ratio  for  the  first  two  modes  are  given  in 
Table  3.1.  These  values  are  very  realistic,  and  the  variance  of  the  damping  is  two  to  four 
times  smaller  than  an  earlier  analysis  made  with  arbitrary,  constant  length  segmentation. 
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Table  3.1  Estimates  of  Natural  Frequency  and  damping  Ratio  for  the  JPL. 


Natural  Frequency  (std.  dev.) 

Damping  Ratio  (std.  dev.) 

First  Mode 

1.02  Hz  (0.007  Hz) 

2.8%  (0.57%) 

Second  Mode 

3.14  Hz  (0.02  Hz) 

3.4  % (0.35%) 

3.3.2  Modeling  of  Earthquake  Strong  Motion 

Earthquakes  can  be  seen  as  high-energy  forcing  functions  that  can  be  used  to  excite 
a system  of  interest.  Since  the  earthquake  is  a one-time  transient  event,  it  is  by  definition 
non-stationary.  It  also  is  an  unique  event,  so  averaging  techniques  available  to  continuous 
excitations  are  not  applicable.  Use  of  parametric  models  have  the  advantage  of  being 
stochastic  in  nature,  parsimonious,  and  directly  related  to  structural  properties.  If  parametric 
models  can  accurately  simulate  a process,  then  the  necessary  filter  has  been  discovered  and 
the  parameters  can  be  used  to  identify  structural  properties.  The  same  reasoning  is  true  if 
the  one-step  predictor  can  be  estimated  with  a parametric  model.  Up  to  this  point,  the 
assumption  has  been  made  that  earthquake  strong  motion  can  be  accurately  modeled  by  AR 
or  ARMA  techniques  without  thoroughly  investigating  whether  this  is  so.  This  important 
assumption  will  now  be  examined. 

In  the  1960’s  and  1970’s  there  was  great  interest  in  modeling  seismic  events  for  the 
purpose  of  identifying  the  source.  The  driving  need  was  the  ability  to  differentiate  nuclear 
explosions  from  seismic  events.  It  was  assumed  that  there  is  a difference  between  the 
ground  motions  due  to  earthquakes  and  explosions.  If  earthquakes  and  explosions  could 
each  be  well  represented  by  an  AR  or  ARMA  model,  the  essence  of  each  process  would  be 
condensed  into  a few  parameters,  greatly  simplifying  the  discrimination  task.  In  addition, 
for  certain  model  structures,  the  parameters  are  related  to  important  physical  properties 
such  as  earthquake  intensity,  duration,  and  distance  (Cakmak  and  Sheriff,  1984). 

A detailed  attempt  to  discriminate  between  blasts  and  earthquakes  derives  an  AR 
model  from  the  initial  displacement  potential  of  the  two  source  mechanisms  (Dargahi- 
Noubary  et  al,  1978).  The  P-wave  history  was  found  to  be  modeled  as  a third-order  AR 
process.  The  parameters  identified  by  this  analysis  are  useful  for  estimating  source 
properties,  such  as  event  duration,  and  reflections  of  the  direct  P-wave.  While  these 
parameters  are  not  of  much  use  to  ascertaining  the  near-surface  soil  properties,  they  are  of 
use  to  the  reflection  seismologist. 

Closer  to  the  problem  at  hand,  a very  understandable  exposition  of  the  segment 
approach  and  the  application  of  the  Burg  estimator  of  AR  parameters  to  earthquake 
characterization  is  given  in  a paper  by  Jurkevics  and  Ulrych  (1978).  In  this  paper  the 
authors  were  interested  in  studying  how  the  process  changed  through  time,  rather  than 
assuming  that  the  system  parameters  were  time-invariant.  The  Orion  Boulevard  recording 
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of  the  1971  San  Fernando  earthquake  was  used  as  the  data  series,  and  the  surface 
accelerogram  modeled  as  a second-order  AR  process  by  the  Burg  method. 

Three  pieces  of  information  are  gleaned  from  the  estimation  — the  two  AR 
parameters  and  the  innovations  variance.  For  the  frequencies  of  interest,  the  system  only 
exhibited  first  mode  response.  The  authors  point  out  that  while  a larger  AR  order  could 
have  been  used  to  model  the  event,  using  the  lowest  order  that  adequately  represents  the 
process  allows  "the  motion  to  be  less  deterministic  in  nature,  since  the  filter  length  is 
shorter."  The  strong  motion  record  is  cut  into  very  short  segments  (one  second,  fifty  data 
points)  since  a short  segment  best  preserves  the  assumption  that  the  process  is  stationary 
within  that  segment.  The  use  of  the  Burg  algorithm  allows  a good  estimation  to  be  made 
from  fifty  points  (see  Chapter  2). 

The  Z-polynomial  representing  the  second-order  system  (Eq.  2.38) 

(3-6) 

is  solved  by  Jurkevics  and  Ulrych  for  the  complex  roots  Zj  and  z^,  which  are  complex 
conjugates  located  just  outside  the  unit  circle.  This  conjugate  pair  is  illustrated  in  Fig.  3.3, 
in  rectangular  coordinates.  Also  shown  are  the  values  given  as  polar  coordinates  R and  0 
(O<0<l/(2At)  Hz).  0 gives  the  frequency  location  of  the  pole,  or  modal  peak,  while  R gives 
the  peakedness  or  half-width  of  the  curve.  The  variance  of  the  innovations  is  the  amplitude 
scaling  factor. 

In  order  to  smooth  the  data,  deterministic  curves  were  fitted  to  the  plotted  discrete 
values  of  the  three  parameters.  The  sharpness  of  the  peak,  R,  changed  somewhat  through 
the  earthquake  process,  while  the  frequency  shift  was  never  more  than  about  eight  percent. 
The  innovation  variance  traced  the  envelope  shape  of  the  accelerogram  through  time.  The 
study  is  quite  successful  in  characterizing  the  time  varying  nature  of  earthquake  strong 
motion.  This  helps  show  the  validity  of  modeling  strong  motion  with  parametric  models. 

An  approach  more  useful  for  identifying  soil  parameters  is  a derivation  of  an  ARMA 
model  for  a SDOF  oscillator  based  on  an  analog  to  the  continuous-time  differential 
equations  of  motion  (Chang  et  al.,  1982).  The  intermediaiy  between  the  continuous  domain 
and  the  discrete  domain  is  the  assumption  of  a match  between  actual  and  calculated 
autocorrelations.  Equations  are  provided  to  calculate  the  damping  ratio  and  resonant 
frequency  of  n-degree-of-freedom  oscillators  from  the  2n  AR  parameters.  Phase  relations 
are  preserved  in  the  MA  parameters. 

Chang  et  al.  used  the  Box-Jenkdns  estimation  procedure  to  calculate  the  ARMA 
parameters  (Box  and  Jenkins,  1976).  The  quality  of  the  estimation  is  decided  by  how  well 
the  residuals  approach  white  noise,  and  is  given  by  the  "Q"  statistic  (Chang  et  al.,  1982). 
The  Box-Jenkins  method  is  only  applicable  to  stationary  time  series,  so  a problem  exists  in 
applying  it  to  earthquake  records.  The  adopted  solution  was  to  divide  the  complete  record 
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Fig.  3.3  Locations  of  filter  poles  Zj  and  in  the  Z plane,  relative  to  the  unit  circle  (from 
Chang  et  al.,  1982). 


57 


into  constant-length  quasi-stationary  segments,  with  five  second  segments  (250  data  points 
each)  giving  the  most  satisfactory  results. 

Strong  motion  records  from  1942  El  Centro,  and  the  1971  San  Fernando  earthquakes 
were  analyzed.  Almost  all  the  data  segments  were  best  modeled  as  ARMA(4,1)  processes. 
The  remaining  segments  were  best  fit  by  an  ARMA(2,1)  model.  For  the  SW  component  of 
the  1940  Imperial  Valley  earthquake,  recorded  at  El  Centro,  and  modeled  as  ARMA(2,1), 
the  parameters  corresponded  to  the  resonant  frequency  ranging  from  2.4  to  3.2  Hz.  The 
corresponding  damping  ratio  ranged  from  27%  to  60%.  However,  it  is  not  clear  to  what 
physical  entity  these  parameters  correspond  to  since  there  is  not  a well-defined  system  — a 
known  input  into  a column  of  soil  (filter)  and  a measured  output.  The  confusion  as  to  what 
physical  quantity  is  being  modeled  is  the  inverse  of  that  discussed  in  Chapter  2.  For  the  case 
under  discussion,  it  appears  that  the  AR  coefficients  model  the  shape  of  the  ground  motion 
history,  and  the  corresponding  spectral  estimate  is  descriptive  of  the  waveform  shape  rather 
than  the  mechanical  behavior  of  the  earth  system. 

The  Box-Jenkins  method  used  by  Chang  et  al.  selects  the  model  order  based  on  the 
data  itself,  rather  than  by  a priori  decision.  The  authors  therefore  expect  that  other 
earthquakes  might  yield  ARMA  orders  different  from  (4,1)  and  (2,1).  The  authors  also 
believe  that  the  source  of  non-stationarity  in  the  surface  strong  motions  is  a result  of  changes 
in  the  input  driving  function  rather  than  in  the  soil  (filter),  a conclusion  with  which  a 
seismologist  would  heartily  agree. 

Segmentation  of  non-stationary  strong  motion  records  is  still  in  use  today.  This 
approach  offers  a straight-forward  way  to  apply  familiar,  tractable  methods  to  a strictly 
inapplicable  problem.  A recent  application  applied  the  AlC-based  segmentation  scheme 
used  by  Gersch  and  Brotherton  (1982)  to  earthquake  records  from  Romania  (Popescu  and 
Demetriu,  1990).  This  choice  was  made  after  a thorough  review  of  the  various  methods 
available,  including  Kalman  filter  methods  that  are  designed  for  non-stationary  data.  An 
ARMA  model  is  fitted  using  the  Maximum  Likelihood  method  with  the  goal  of 
characterization  and  simulation  of  earthquakes  with  few  parameters. 

A preliminary  canonical  correlation  analysis  provides  an  "good"  initial  guess  as  to  the 
ARMA  orders,  and  the  final  choice  is  made  from  statistics  developed  during  the  Maximum 
Likelihood  determination  and  application  of  the  AIC.  For  simulation,  each  segment  can  be 
calculated  by  passing  white  noise  with  variance  equal  to  that  of  the  associated  ARMA  model 
through  the  calculated  process  filter,  and  then  multiplying  the  section  with  an  "envelope" 
function  taken  from  the  actual  intensity  of  the  associated  earthquake  segment. 

In  general,  the  simulations  of  Popescu  and  Demetriu  matched  the  actual  motions  very 
well,  with  statistics  such  as  cumulative  energy  and  root  mean  square  acceleration  tracking 
throughout  the  record.  However,  there  were  some  problems  with  the  match  for  short  term 
energy  for  two  directional  components.  Short-term  spectra  were  also  calculated  for  each 
quasi-stationary  segment,  giving  a graphical  image  of  the  change  in  the  shaking  process 
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through  time.  This  paper  illustrates  the  utility  of  the  segmentation  method  of  approaching 
non-stationary  signals,  and  shows  that  earthquake  strong  motion  can  be  well  modeled  by  the 
ARMA  scheme. 


3.3.3  Transforming  Non-stationary  Signals  into  Stationary  Signals 

Another  approach  to  analyzing  non-stationary  data  using  methods  strictly  stationary 
in  derivation,  is  to  transform  the  non-stationary  data  series  into  a stationary  record.  The 
oldest  and  simplest  technique  is  to  apply  the  difference  operator  to  the  data.  Often  times 
the  new  data  series  will  be  stationary  and  the  parameters  solved  for  with  customary 
techniques.  The  parameters  can  be  adjusted  to  take  into  account  the  effect  of  the 
differencing  transformation  (Bendat  and  Piersol,  1986). 

A more  advanced  approach  is  based  on  the  same  assumption  made  by  Gersch  and 
Brotherton  (1982),  that  the  non-stationary  signal  is  a linear  combination  of  two  components; 

(1)  a stationary  component  that  is  a result  of  the  response  of  the  structure  in  question,  since 
the  physical  properties  are  generally  thought  to  be  time  invariant  (a  problem  in  soils  for 
strains  greater  than  the  threshold  strain),  and  (2)  the  second  component  is  a non-stationary 
auxiliary  input  to  time  varying  noise  processes. 

Based  on  the  loose  definition  of  stationarity  used  in  this  report  (data  variance  is 
constant  through  time),  this  approach  sees  the  earthquake  strong  motion  as  a constant-mean 
process  acted  on  by  an  "envelope"  function  having  time-varying  amplitude.  Polhemus  and 
Cakmak  (1981)  construct  a polynomial  to  describe  the  change  in  variance  over  time  and 
employ  it  to  "correct"  the  variance  to  a constant.  They  use  this  method  to  simulate 
earthquake  strong-motion  using  an  ARMA  model  (1981).  A particular  earthquake  is  seen 
as  an  instance  of  a stochastic  process.  Therefore,  a "similar"  earthquake  would  not  have  to 
be  identical,  but  only  need  to  have  the  same  stochastic  descriptors.  Indeed,  this  is  the 
stochastic  foundation  that  makes  the  AR  and  ARMA  methods  possible. 

Polhemus  and  Cakmak  (1981)  devised  a three-stage  procedure  to  estimate  constant 
parameters  from  a time  series  with  a time-varying  variance: 

(1)  An  estimate  of  the  variance  function  a/(t)  is  made  in  a non-biased  manner.  A "power 

transformation"  is  used  and  fitted  to  a polynomial  using  least  squares. 

(2)  The  acceleration  series  is  converted  to  a stationary  series.  The  standard  deviation  at 

each  time  step  is  derived  from  step  (1)  and  is  used  to  normalize  the  non-stationary 

accelerations  history  into  a constant-variance  data  series. 

(3)  An  ARMA  model  is  fit  to  the  transformed  data  and  the  constant  process  parameters 

estimated. 
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A new,  statistically  equivalent,  earthquake  record  can  be  simulated  by  passing  white 
noise  through  the  derived  ARMA  model  and  then  "denormalizing"  the  signal  with  the  step- 
by-step  standard  deviations  calculated  in  step  (2).  The  model  is  validated  by  examination 
of  the  residuals  of  the  actual  and  simulated  data.  If  the  residuals  from  the  two  are 
uncorrelated  (white),  then  the  model  has  accurately  represented  the  process. 

This  method  has  been  applied  to  a suite  of  earthquake  records  in  an  attempt  to 
estimate  seismological  variables  (Cakmak  and  Sherif,  1984).  The  model  parameters  are  used 
to  estimate  values  for  temblor  intensity,  duration,  distance  from  the  site  to  the  fault,  and 
local  geology.  Statistical  methods  are  employed  to  fit  these  variables  to  model  parameters, 
which  now  include  parameters  describing  the  variance  function  a^(i)  as  well  as  the  familiar 
ARMA  parameters.  The  brief  results  for  the  1971  San  Fernando  earthquake  are  positive, 
but  the  authors  note  that  much  more  work  needs  to  be  done  before  the  assumptions  can  be 
accepted. 


3.3.4  Modeling  of  Non-stationary  Processes  — Time  Adaptive  filtering 

A non-stationary  process  can  be  modeled  directly  using  time  adaptive  filtering,  where 
the  estimates  of  the  AR  parameters  are  updated  for  each  time  step.  An  early 
implementation  of  this  approach  to  earthquake  strong  motion  data  was  in  the  paper  by 
Jurkevics  and  Ulrych  (1978)  discussed  above.  The  authors  modified  the  Least  Mean  Squares 
algorithm  (Alexander,  1986)  to  operate  on  the  data  forward  and  backward,  after  Burg 
(1975).  This  modification  results  in  instantaneous  spectral  estimates  with  shapes  similar  to 
that  calculated  by  applying  the  Burg  method  to  segmented  data  (see  Sect.  3.3.2).  In 
addition,  the  step  size  of  the  adaptive  algorithm  changed  as  a function  of  the  local  input 
amplitude. 

The  modified  algorithm  of  Jurkevics  and  Ulrych  allows  the  filter  to  be  applied  to  the 
data  in  reverse,  in  order  to  eliminate  startup  error.  Startup  error  due  to  poor  initial 
estimates  and  a small  past  history  is  especially  difficult  with  seismic  signals  which  are 
minimum  phase,  i.e.  energy  density  is  concentrated  at  the  beginning  of  the  signal.  The 
instantaneous  estimates  of  the  two  AR  parameters  and  the  innovations  variance  were 
averaged  in  fifty  data  point  segments  to  smooth  the  noise,  and  deterministic  curves  fitted  to 
the  data.  The  proper  time  step  to  use  in  the  adaptive  algorithm  was  somewhat  troublesome 
to  determine,  an  improper  value  either  critically  underdamps  or  overdamps  the  response. 

The  adaptive  estimations  reported  by  Jurkevics  and  Ulrych  were  virtually  identical 
to  those  calculated  by  the  segmentation  method.  This  result  shows  that  carefully  segmenting 
non-stationary  data  does  a very  good  job  of  characterizing  the  local,  quasi-stationary  process. 
The  frequency  and  amplitude  content  of  the  accelerograms  simulated  by  the  resulting  filters 
were  the  same  as  the  actual  recorded  signal. 
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A much  more  direct  approach  is  to  utilize  the  Kalman  filter  estimation  of  the  AR 
parameters.  Kitagawa  extended  the  work  of  Akaike  through  a state-space  AR  model  with 
constrained  time-varying  coefficients  (1983).  The  model  assumes  a zero  mean,  non- 
stationary data  series,  such  as  an  accelerogram.  In  contrast  to  the  method  of  Jurkevics  and 
Ulrych  discussed  above,  where  the  AR  coefficients  could  change  abruptly  through  time  and 
then  had  to  be  smoothed  and  filtered  by  an  ad  hoc  technique,  Kitagawa  employs  a linear 
constraint  to  insure  that  the  AR  coefficients  can  only  change  gradually  through  time. 

The  time  varying  AR  model  is  defined  as 

m 

+ (3-7) 

j=i 

where  m = counter  for  AR  order  l..k 

n = time  step  1..N  number  of  data  points. 

where  the  AR  coefficient  is  a gradual  function  of  the  time  step  counter  n.  A simple  least 
squares  fit  would  entail  the  calculation  of  mN  parameters  from  N data  points.  This  large 
overparameterization  leads  to  poor  estimates.  A parsimonious  approach  is  to  treat  the  AR 
coefficients  as  random  variables  and  apply  a stochastic  rule  as  to  their  behavior  through 
time.  The  change  in  the  AR  parameters  is  constrained  by  a "stochastically  perturbed 
difference  equation",  which  must  reflect  the  actual  behavior  of  the  process  being  modeled. 
The  estimation  is  now  a constrained  least  squares  problem. 

For  the  simplest  case,  the  parameters  are  modeled  as  randomly  changing,  modeled 
by  a random  walk  prior  distribution,  and  neighboring  coefficients  are  almost  equal.  The 
equation  governing  the  change  is  a first  order  difference  model 

a(i,n)=a(i^n-l)+6(i,n) 

where  6 = Gaussian  white  noise  with  0 mean,  variance 

n = data  point  number  (time  step). 

If  the  process,  hence  AR  parameters,  changes  more  slowly,  the  constraint  is  better  defined 
by  a second  order  difference  model 

a(i,  n)  = 2a(i,n-l)-a  (z,  n -2) +6  (z,  n). 

The  order,  k,  of  the  difference  equation  that  best  models  how  rapidly  the  system  is  changing 
through  time,  must  be  estimated,  and  are  parameters  of  the  model  along  with  the  AR  order 
m and  innovations  variance.  The  variance  of  the  constraint  innovations,  controls  the 
"intensity"  of  the  stochastic  change,  and  must  also  be  estimated. 
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There  are  now  (m+k+2)  parameters  to  be  estimated  — k AR  parameters,  m 
constraints,  and  the  two  variances,  and  r^,  which  can  be  combined  into  the  ratio  /i  = 
The  ratio  ii  represents  the  trade-off  between  the  (in)fidelity  to  the  actual  data  and 
the  (in)fidelity  to  the  change  constraint.  The  trade-off  ratio  has  been  likened  to  a signal-to- 
noise  ratio,  or  in  Bayesian  terms  a measure  of  belief  in  the  a priori  assumptions.  Non-linear 
maximization  of  the  likelihood  parameter  is  used  to  estimate  /i,  and  the  AIC  is  used  to 
define  the  system  parameters  m,  and  k.  At  this  point  Kalman  filtering  is  used  to  estimate 
the  values  of  the  model  parameters. 

An  instantaneous  power  spectral  density  can  be  generated  from  the  "instantaneous" 
AR  estimates.  Given  the  local  AR  parameters  just  estimated,  a(j,n)  and  the  spectral 
estimate  can  be  defined  by 


5(o),«)  = 


o^(n) 


;=i 


(3.10) 


where  i 
m 
n 

G) 

j 


= y-i 

= AR  order 
= time  step  of  interest 
= radial  frequency 
= counter  for  l..m  AR  order 


This  spectrum  is  simply  a helpful  visual  aid  as  to  the  behavior  of  the  process  at  that  point 
in  time. 


This  entire  scheme  was  tested  by  comparing  the  theoretical  and  estimated 
instantaneous  spectra  for  synthesized  non-stationary  AR  processes.  An  example  slowly- 
changing  time  series  is  shown  in  Fig.  3.4.  The  actual  gradually  changing  spectrum  for  this 
data  is  shown  in  Fig.  3.5,  and  should  be  compared  to  the  estimated  spectrum  shown  in  Fig. 
3.6.  The  estimated  spectra,  for  AR  order  5,  second  order  difference  constraint,  and  trade- 
off ratio  9.6x10'^,  is  virtually  the  same  as  the  theoretical  solution.  The  match  for  rapidly  and 
instantaneously  changing  processes  match  almost  as  well. 

An  earlier  paper  by  Kozin  (1977)  presented  a very  similar  approach.  However,  Kozin 
used  an  orthogonal  family  of  Legendre  polynomials  to  "constrain"  the  changing  AR 
coefficients.  He  also  used  a consistency  method  to  apply  maximum  likelihood  to  non- 
stationary  signals  and  AIC  to  determine  the  best  choice  for  model  order. 

The  use  of  Kalman  filtering  is  presented  in  more  detail  in  later  papers,  where  the 
method  is  directly  applied  to  characterizing  earthquake  strong  motion  records  (Kitagawa  and 
Gersch,  1985;  Gersch  and  Kitagawa,  1985).  The  authors  point  out  that  earthquake 
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acceleration  history  is  modulated  by  a slowly  changing  envelope  function,  as  described  by 
Polhemus  and  Cakmak  above.  This  envelope  is  the  change  in  scaling  of  the  AR  parameters 
through  time,  or  the  "smoothed  value  of  the  instantaneous  variance."  The  Kalman-Akaike 
based  technique  just  discussed  (Kitagawa,  1983)  can  be  used  to  estimate  the  smoothed  trend 
of  the  innovations  variance,  or  changing  power. 

The  "smoothness  priors  time-varying  AR  coefficient"  model  was  applied  to  the  1971 
San  Fernando  earthquake  (Gersch  and  Kitagawa,  1985)  and  the  1949  Olympia,  Washington 
earthquake  (Kitagawa  and  Gersch,  1985).  The  results  for  the  Olympia  temblor  will  be  used 
as  an  example.  The  authors  recalculated  and  optimized  model  and  constraint  order  every 
one-hundred  data  points,  out  of  a total  twelve-hundred.  The  actual  accelerogram  is  shown 
in  Fig.  3.7,  with  a simulated  history  shown  in  Fig.  3.8  for  comparison. 

The  change  of  the  system  through  time  is  given  by  the  change  in  AR  coefficients 
through  time,  shown  in  Fig.  3,9,  and  the  innovations  variance  history  shown  in  Fig.  3.10. 
Notice  that  the  system  descriptors  change  quite  infrequently,  and  with  "quantum"  jumps. 
The  variance  envelope,  or  power,  shows  several  "bumps"  which  the  authors  believe  are 
indicative  of  the  arrival  of  the  various  wave  modes.  The  local  impulse  response  function  of 
the  process  can  be  computed  from  the  instantaneous  AR  coefficients  and  innovations 
variance.  This  information  is  combined  in  Fig.  3.11a,  which  shows  the  instantaneous  spectral 
description  of  the  system  through  time.  The  spectra  actually  reveal  no  more  information 
than  the  other  plots,  but  presents  the  physical  relationships  in  a more  familiar  manner.  For 
comparison.  Fig.  3,11b  shows  the  instantaneous  spectrum  for  the  simulated  time  history, 
shown  in  Fig.  3.8, 

At  this  writing,  these  methods  have  only  been  applied  to  model  the  strong  motion 
response.  It  is  hoped  that  Gersch  and  Kitagawa  will  adapt  this  method  to  estimating  the 
parameters  of  a single  input-single  output  system.  Attempts  are  currently  being  made  to 
apply  extended  Kalman  filters  to  directly  identifying  nonlinear  soil  properties  (Lin,  1990). 
This  preliminary  study  was  made  for  a simulated  soil  system  for  which  the  nonlinearity  was 
assumed  to  be  purely  hysteretic.  The  input  and  output  signals  for  the  soil  system  are 
required  for  parameter  identification,  and  estimates  of  the  statistical  nature  of  the  noise. 

The  results  were  varied,  with  fair  success  when  the  model-generated  backbone  curves 
were  the  desired  output.  If  the  soil  system  is  well  characterized  by  the  chosen  model,  then 
only  one  input-output  data  series  is  needed  to  completely  describe  the  soil’s  nonlinear 
behavior.  If  not,  strong  motion  records  for  several  different  strain  levels  are  needed.  At  this 
point  many  assumptions  were  made  and  the  system  limited,  but  the  results  were  heartening 
and  further  work  on  these  lines  is  warranted. 
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Fig.  3.4  An  artificially  generated  AR  process  with  a gradually  changing  spectrum 
(from  Kitagawa,  1983) 
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Fig.  3.5  Actual  instantaneous  spectra  for 
time  series  shown  in  Fig.  3.4 
(from  Kitagawa,  1983). 
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Fig.  3.6  Calculated  instantaneous  spectra 
for  time  series  shown  in  Fig.  3.4 
(from  Kitagawa,  1983). 


Fig.  3.7  Actual  accelerogram  from  the  1949  Olympia,  Washington  earthquake  (from 
Kitagawa  and  Gersch,  1985). 


Fig.  3.8  Simulated  time  history  for  the  1949  Olympia,  Washington  earthquake  (from 
Kitagawa  and  Gersch,  1985). 
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Fig.  3.9  Time  varying  AR  model  coefficients  (A)  a(l,l),  (B)  a(2,l),  (C)  a(3,l)  for  the  1949 
Olympia  temblor  (from  Kitagawa  and  Gersch,  1985). 


Fig.  3.10  Smoothed  instantaneous  innovations  variance  for  the  1949  Olympia  temblor  (from 
Kitagawa  and  Gersch,  1985). 


Fig.  3.11a  Actual  instantaneous  spectrum 
for  1949  Olympia  earthquake  (from 
Kitagawa  and  Gersch,  1985). 


Fig.  3.11b  Instantaneous  spectrum  for  the 
simulated  1949  Olympia 
earthquake  (from  Kitagawa  and 
Gersch,  1985). 
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3.3.5  Summary 


Several  different  methods  have  been  presented  to  directly  assess  non-stationary  data. 
These  methods  cover  a wide  range,  from  transforming  the  data  into  a stationary  signal  to 
full-fledged  non-linear,  non-stationary  analysis.  The  method  used  will  depend  on  the  nature 
of  the  data  available  and  the  nature  of  the  information  needed  to  be  gleaned.  Transforming 
the  data  (Polhemus  and  Cakmak,  1981)  is  an  acceptable  method  if  the  non-stationarity  is 
mild  and  mostly  in  the  amplitude.  This  method  does  not  give  much  information  as  to  the 
time-varying  change  in  frequency  and  damping. 

Segmentation  of  the  data  into  stationary  pieces  gives  very  good  results  and  is  widely 
used.  The  approach  by  Gersch  and  Brotherton  (1982)  is  the  most  rational.  It  makes  the 
best  use  of  the  limited  data  available  and  allows  the  noise  model  to  change  through  time. 
However,  if  the  process  of  interest  has  time  varying  parameters,  then  a recursive  technique, 
such  as  presented  by  Kitagawa  (1983),  must  be  used.  In  fact,  the  extended  Kalman  filter  is 
often  used  in  other  fields  to  estimate  smoothly  varying  non-linear  parameters. 


3.4  Estimation  of  Soil  Parameters  From  Earthquake  Strong  Motion  Data 
3.4.1  Introduction 

Up  to  this  point,  the  work  presented  has  been  of  a theoretical  bent,  to  show  the 
methods  available  for  estimating  system  parameters.  Understanding  these  discussions  will 
allow  the  reader  to  understand  actual  field  applications.  In  this  section,  actual  field 
experiments  made  using  earthquake  strong  motion  data  are  presented. 


3.4.2  Use  of  Shear  Beam  Theory 

An  early  success  using  measured  response  of  an  earth  dam  to  earthquake  excitation 
to  estimate  soil  properties  was  made  for  the  Santa  Felicia  dam  in  California  (Abdel-Ghaffar 
and  Scott,  1979).  The  input  (base  displacement)  and  output  (crest  motion)  of  the  structure 
was  recorded  for  the  1971  San  Fernando  earthquake,  and  another  temblor  in  1976  (Ml  = 
3.7).  A shear-beam  theory  model  was  used  to  estimate  shear  moduli  and  damping  ratio  for 
the  rolled-fill  earth  dam,  constructed  from  gravelly  sands. 

The  soil  was  modeled  as  a hysteretic  SDOF  system  with  a non-linear  restoring  force. 
An  estimate  of  the  hysteretic  response  can  be  made  by  plotting  the  crest  displacement 
relative  to  the  base,  against  absolute  acceleration  of  the  dam.  Using  the  shear-beam  model, 
the  slope  and  area  of  the  hysteresis  loop  can  correspond  to  shear  modulus  and  damping 
ratio  (see  Richart  et  al,  1970).  Correlations  can  also  be  made  to  a finite  element  model 
using  estimated  parameters.  The  fundamental  frequency  behavior  of  the  dam  was  isolated 
by  narrow  band-pass  filtering  the  data.  Without  this  filtering,  the  plotted  hysteresis  loops 
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were  erratic  and  useless  for  analysis.  The  precariousness  of  this  modeling  can  be 
appreciated  from  the  very  small  displacement  of  the  top  of  the  dam  relative  to  the  base  — 
2.0  mm. 

The  non-stationary  response  of  the  dam  was  analyzed  by  breaking  the  record  into 
approximately  one-second  segments.  Using  both  two-dimensional  shear-beam  theory  (with 
constant  shear  modulus)  and  one-dimensional  theory  (shear  modulus  varying  with  depth), 
a relationship  for  dynamic  stiffness,  G,  was  derived.  For  the  two-dimensional  case  the 
relationship  for  G is  given  by  : 


. (3.11) 


where 

= fundamental  mode  shear  stress 

Yn 

= fundamental  mode  shear  strain 

P 

= mass  density  of  dam  soil 

Vs 

= shear  velocity  of  dam  soil 

(Oil 

= fundamental  mode  natural  frequency 

Y 

-^max 

= maximum  relative  displacement  of  hysteresis  loop. 

= maximum  absolute  acceleration  of  hysteresis  loop 

A similar  relationship  exists  for  the  one-dimensional  case.  In  addition,  expressions  for  stress 
and  strain  at  any  time  were  derived. 

With  these  equations  the  relationship  between  shear  modulus  and  shear  strain  can 
be  calculated.  Additional  information  is  given  by  the  data  from  the  second,  smaller 
earthquake.  An  example  of  the  results  is  shown  in  Fig.  3.12,  where  the  stiffness  reduction 
values  calculated  are  well  within  the  range  of  those  expected.  However,  a similar  analysis 
for  the  damping  ratio  yielded  a curve  of  different  shape  than  expected  — S-shaped  rather 
than  hyperbolic,  as  shown  in  Fig.  3.13.  The  discrepancy  might  be  due  to  the  fact  that  the 
soils  tested  were  predominantly  gravels,  and  there  have  been  virtually  no  laboratory  tests  run 
on  gravel,  and  the  results  were  accurate.  However,  it  has  been  seen  throughout  Chapt.  2 
and  3 that  the  estimations  for  damping  are  often  inaccurate  and  inexact. 
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Fig.  3.12  Experimentally  measured  modulus  degradation  compared  to  previous  values  (from 
Abdel-Ghaffar  and  Scott,  1979). 


Fig.  3.13  Experimentally  measured  damping  ratios  (from  Abdel-Ghaffar  and  Scott,  1979). 
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3.4.3  Fourier  Analysis  of  Non-stationary  Ground  Motion 

A study  has  been  made  of  data  recorded  at  the  Wildlife  site  in  the  Imperial  Valley 
during  the  1987  Superstition  Hills  earthquake  (Vucetic  and  Zorapapel,  1990).  For  this 
event,  three-dimensional  accelerograms  were  recorded  at  a depth  of  7.5  meters  at  the 
surface.  In  addition,  pore  water  pressure  was  monitored  throughout  the  shaking.  The  paper 
attempts  to  correlate  a decrease  in  natural  frequency  of  the  soil  layer,  associated  with 
degradation  of  the  shear  modulus,  with  the  increase  of  pore  water  pressure  as  the  condition 
of  liquefaction  is  approached.  However,  it  has  been  postulated  that  the  piezometers  were 
not  functioning  properly  (Hushmand  et  al.,  1992,  1991).  The  methods  used  in  the  analysis 
are  instructive,  in  addition  to  giving  a view  of  the  state-of-the-practice  in  the  use  of  system 
identification  in  the  geotechnical  field. 

The  system  is  well-defined  for  this  problem,  since  both  the  input  and  output  are 
known.  Rather  than  use  the  more  certain  AR  methods  discussed  in  Chapt.  2,  the  authors 
choose  to  characterize  the  motions  by  simple  Fourier  spectra,  with  all  the  attendant 
uncertainty  and  problems  discussed  in  the  previous  chapter.  The  system  response  function 
is  not  calculated  directly  through  a least  squares  deconvolution  or  system  identification 
procedure,  rather  is  characterized  by  the  spectral  ratio.  Because  it  is  a simple  ratio,  a small 
error  in  one  of  the  spectral  estimates  can  have  a very  large  effect  on  the  calculated  value. 

The  non-stationarity  of  the  data  is  taken  into  account  by  segmenting  the  data, 
although  the  segmenting  is  not  based  on  preserving  sections  with  a common  variance. 
Instead,  the  segmentation  is  based  on  physical  concerns  such  as  arrivals  of  various  wave 
modes,  or  changes  in  pore  water  pressure.  The  analysis  showed  that  the  fundamental  period 
of  the  soil  layer  lengthened  as  the  pore  pressure  built  up  and  the  soil  stiffness  degraded. 
The  entire  process  was  also  modeled  using  the  DESRAMOD  computer  program  (Vucetic, 
1986)  which  calculates  soil  and  pore  water  behavior  based  on  a one-dimensional  lumped 
parameter  model.  The  resonant  frequencies  calculated  by  DESRAMOD  and  those  from  the 
Fourier  analysis  agreed  well,  especially  for  the  vertical  component. 


3.4.4  Parameter  Estimation  Using  Impedance  Functions 

Actual  earthquake  excitation  has  been  used  to  measure  experimental  impedance 
functions  of  structures  for  comparison  to  theoretical  functions,  with  acceptable  results  (Mau 
and  Wang,  1990;  Loh  and  Mau,  1989).  Loh  and  Mau  studied  the  modal  behavior  of  the 
one-quarter  scale  model  nuclear  power  plant  containment  vessel  in  Lotung,  Taiwan.  The 
rigid  structure  is  very  well  instrumented  so  that  the  rigid  body  motion  could  be  very  well 
characterized  with  simple  computation.  The  paper  reports  preliminary  work,  for  which  the 
free-field  motion  is  not  taken  into  account.  The  best  fit  between  the  theoretical  impedance 
functions  and  the  experimental  impedance  functions  was  for  the  assumed  material  values 
given  in  Table  3.2. 
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Table  3.2  Estimated  Material  Parameters  for  Lotung  Site  (after  Loh  and  Man,  1989), 


Natural  Frequency 

2.72  Hz 

Damping  Ratio 

10.9  % 

Shear-wave  Velocity 

310  m/s 

Man  and  Wang  (1990)  continue  the  work  on  the  Lotung  containment  structure  with 
a more  complete  analysis  where  the  form  of  the  impedance  function  was  not  assumed. 
Again,  only  the  rocking  impedance  was  studied  to  simplify  the  embedded  foundation 
problem.  The  coupling  terms  are  ignored,  allowing  the  system  to  be  modeled  as  a single- 
input-single-output system,  with  impedance  function  acting  as  the  transfer  function  (filter). 
The  kinematic  input  motions  must  then  be  calculated  from  the  measured  motions. 

The  dynamic  impedance  is  calculated  using  the  method  of  Bendat  and  Piersol  (1986), 
discussed  in  Chapt.  2,  Sect,  2.2.1.  The  transfer  function  is  calculated  from  the  ratio  of  the 
cross-spectrum  to  the  autospectrum  (Eq.  2.3).  The  coherence  function  was  also  calculated, 
and  only  relationships  for  frequencies  with  high  coherence  (1.5  - 4.5  Hz)  were  used  to  define 
the  impedance  function.  The  results  were  very  poor,  possibly  because  the  assumptions  made 
about  lack  of  coupling  were  incorrect,  or  more  than  one  mode  must  be  included  in  the 
analysis  to  obtain  useful  results. 


3.4.5  Non-linear  Ground  Response  Analyses 

Possibly  the  best  set  of  data  for  earthquake  excitation  of  soils  exists  for  the  test  site 
operated  by  the  Electric  Power  Research  Institute  (EPRI)  and  the  Taiwan  Power  Co.  at 
Lotung  Taiwan  (Tang,  1987).  At  this  site,  two  locations  are  instrumented  with  three- 
component  accelerometers  at  depths  of  47,  17,  11,  6 meters,  and  at  the  surface.  One  array 
is  in  the  free-field  while  the  other  is  adjacent  to  the  one-quarter  scale  nuclear  containment 
vessel  mentioned  in  Section  3.4.4.  The  site  is  also  well  instrumented  with  piezometers  at 
various  depths  and  locations.  The  simplified  soil  profile  consists  of  30-35  m of  silty  sand  and 
sandy  silt  with  some  gravel,  above  clayey  silt  and  silty  clay.  The  water  table  is  within  half 
a meter  of  the  ground  surface.  This  area  is  seismically  active,  and  many  earthquakes  of 
exhibiting  a wide  range  of  magnitudes  have  occurred  since  1986. 

A series  of  studies  have  been  undertaken  at  Lotung  using  parameter  identification 
to  evaluate  the  non-linear  response  of  soils  due  to  strong  motion  (Chang  et  al.,  1991,  1990, 
1989).  The  authors  note  that  outside  of  actual  liquefaction  sites,  there  has  been  no 
documented  field  evidence  of  the  degradation  of  soil  properties  with  increasing  strain  (Chang 
et  al,,  1989).  One  notable  exception  is  the  work  done  for  the  U.S,  Nuclear  Regulatory 
Commission  (Shannon-Wilson  and  Agbabian,  1975). 
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The  first  phase  of  the  work  concentrated  on  calculating  a transfer  function  for  the  soil 
at  several  different  points  of  the  excitation  history.  The  parameter  of  interest  from  the 
transfer  function  was  resonant  frequency.  A simphfied  analysis  for  a two-layer  system  was 
used  to  calculate  the  effective  S-wave  velocity  since  the  soil  profile  was  known  (Dobry  at  al., 
1976).  Shear  moduli  were  calculated  from  the  estimates  of  wave  velocity.  The  strong 
motion  record  was  divided  in  three  sections:  (1)  initial  motion  before  strain  levels  high 
enough  to  degrade  the  soil  was  reached,  (2)  peak  motion,  and  (3)  the  coda  or  ring-down. 

The  transfer  functions  were  estimated  from  the  ratios  of  output  Fourier  spectra  and 
input  spectra  fi-om  various  depths.  The  Fourier  spectra  were  smoothed  by  a triangular  lag 
window.  Unfortunately,  the  data  lengths  were  very  short  since  two  of  the  time  windows  were 
only  four  seconds  long.  A Burg  estimator  would  be  much  more  appropriate  for  this 
application.  In  addition,  use  of  an  optimal  segmentation  scheme  (Gersch  and  Brotherton, 
1982)  might  be  very  helpful  in  identi^ng  exactly  when  the  parameters  of  the  system  actually 
changed. 

The  shear  modulus  showed  a marked  decrease  from  the  initial  low-level  excitation 
to  peak  deformation.  Depending  on  the  depth  used  as  input,  the  modulus  reduction  factor 
ranged  from  0.60  to  0.14.  These  values  are  for  a magnitude  6.5  earthquake  exhibiting  a 
peak  horizontal  acceleration  of  0.21  g.  As  a check  on  reality,  the  S-wave  amplitudes  for  the 
initial  segment  were  in  very  good  agreement  with  the  actual  measured  values.  A 
troublesome  point  is  the  large  discrepancy  in  modulus  reduction  factor  for  the  two  horizontal 
components.  The  expectation  would  be  that  they  would  be  virtually  identical,  since  soil 
degradation  is  usually  thought  of  as  a scalar  quantity. 

This  body  of  data  is  ideal  for  checking  the  results  of  ground  response  programs  such 
as  SHAKE  (Schnabel  et  al,  1972)  and  DESRA-2  (Lee  and  Finn,  1985).  SHAKE  was  tested 
in  both  the  forward  prediction  mode  and  the  inverse  surface-to-depth  mode,  and  DESRA-2 
in  the  forward  mode  (Chang  et  al.,  1990).  The  modulus  degradation  curve  was  calculated 
from  actual  field  data  in  the  manner  discussed  immediately  above.  The  equivalent  damping 
ratio  was  estimated  from  resonant  column  tests  and  the  Seed-Idriss  curve.  For  the  SHAKE 
analyses,  the  additional  input  parameter  was  the  field  gathered  S-wave  velocity  profile. 

The  results  from  SHAKE  show  that  the  calculated  motions  are  higher  than  actual  for 
frequencies  greater  than  0.6  — 1.5  Hz,  with  the  lower  frequency  associated  with  the  analysis 
of  a thicker  layer  of  soil.  There  was  no  correlation  for  phase  information.  The  backward 
analysis  yielded  better  results  than  the  forward  propagation  analyses.  The  disagreements 
might  be  due  to  the  equivalent  linear  use  of  a single  value  of  shear  modulus  per  layer.  The 
modulus  used  might  only  be  valid  for  part  of  the  load  history,  most  specifically  the  peak 
strains  for  which  the  modulus  is  corrected  for  (Chang  et  al.,  1990). 

The  results  from  the  forward  propagation  non-linear  DESRA-2  analysis  show  good 
agreement  between  actual  and  calculated  displacement  for  frequencies  up  to  about  6 Hz. 
There  was  also  good  correlation  for  phase  angle  for  frequencies  up  to  about  3 Hz.  The 
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shape  of  the  equivalent  damping  ratio  curve  back-calculated  from  the  non-linear  analysis  is 
different  from  that  normally  expected.  Rather  than  the  expected  hyperbolic  curve,  the  field 
curve  is  S-shaped.  This  uncommon  S-shaped  curve  is  the  same  shape  as  that  calculated  by 
Abdel-Ghaffar  and  Scott  (1979)  from  actual  field  data  (see  Section  3.4.2). 

Further  work  on  the  shear  modulus  reduction  curve  based  on  actual  large  strain 
measurement  is  reported  by  Chang  et  al.  (1991).  This  study  utilized  seven  different 
earthquakes  with  magnitudes  ranging  from  4.5  to  7.0,  and  peak  horizontal  accelerations 
ranging  from  0.03  g to  0.21  g.  As  for  the  earlier  report  discussed  above  (Chang  et  al.,  1989), 
shear  modulus  is  estimated  from  S-wave  velocities  derived  from  resonant  frequencies  and 
soil  layer  thickness.  Shearing  strains  are  calculated  from  wave  propagation  theory  using 
SHAKE. 

In  this  study,  the  non-stationarity  of  the  accelerograms  are  not  taken  into  account, 
with  the  transfer  functions  being  estimated  for  the  entire  time  history  at  once.  This  means 
that  the  resultant  resonant  frequency  estimated  is  an  average  value  for  the  entire  strain 
history  of  the  earthquake.  The  earlier  work  showed  that  there  can  be  a very  substantial 
change  in  soil  stiffness  through  time,  so  this  approach  seems  incomplete.  Of  course,  for 
small  events  with  little  or  no  degradation,  the  seismic  velocities  and  the  back  calculated 
values  match  very  well. 

The  equivalent  stiffnesses  or  S-wave  velocities  calculated  from  the  transfer  functions, 
are  input  into  the  computer  program  SHAKE,  along  with  damping  curves  from  laboratory 
tests.  Shearing  strains  are  calculated  when  recorded  surface  motions  are  run  through 
SHAKE.  The  calculated  strains  are  open  to  doubt  since  the  program  makes  its  estimation 
from  the  peak  stress  value,  and  calculated  strain  and  stress  at  different  points  in  the  soil 
layer.  All  the  problems  with  the  equivalent-linear  approach  brought  up  by  Chang  et  al. 
(1990)  are  appropriate  here  too. 

If  the  described  method  is  accepted  as  valid,  the  resulting  modulus  reduction  curves 
are  as  shown  in  Fig.  3.14,  along  with  a comparison  to  laboratory  test  results.  The  results 
show  that  for  small  strains  (surface  acceleration  < 0.03  g),  the  results  from  geophysical 
methods,  resonant  column,  and  system  identification  are  in  good  agreement.  For 
intermediate  strains  (5x10'^  to  2x10'^)  the  back  calculated  values  for  modulus  reduction  are 
up  to  twenty  percent  lower  than  the  resonant  column  measurements.  For  higher  strains  of 
3x10'^  to  2xl0'\  the  results  from  cyclic  triaxial  tests  are  in  fair  agreement  with  the  field 
values  for  shear  modulus,  with  moderated  scatter  for  the  laboratory  data. 

This  field  test  is  the  most  complete  to  date.  However,  some  compromising 
assumptions  should  be  mentioned.  The  values  of  S-wave  velocity  and  strain  are  average 
values  for  the  entire  event,  from  low  stress  to  peak  and  back,  and  not  actual  values  for  any 
particular  part  of  the  excitation  history.  The  method  of  calculating  the  transfer  function  is 
can  amplify  inaccuracies,  since  the  high  variance  in  Fourier  estimates  can  cause  even  larger 
errors  when  the  ratio  is  taken.  The 
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Normalized  Shear  Modulus,  G/Gmax 


Fig.  3.14  Back-calculated  normalized  shear  modulus  compared  with  laboratory  test  data 
(from  Chang  et  al.,  1991). 
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use  of  SHAKE  to  calculate  strain  raises  questions,  especially  given  the  problems  discussed 
by  the  authors  one  year  earlier  (Chang  et  al.,  1990).  It  would  be  very  informative  to  use  this 
fine  set  of  data  with  a more  thorough  and  challenging  set  of  analyses. 

Due  to  the  general  lack  of  data  sets  for  which  the  ground  motion  input  into  a soil 
profile  is  known  from  downhole  instruments,  an  attempt  has  been  made  to  simplify  the 
inverse  problem  so  that  only  the  output  (surface)  strong  motion  needs  to  be  known 
(Tokimatsu  et  al.,  1989).  As  in  the  attempt  by  Nakamura  (1989),  an  assumption  has  to  be 
made  as  to  the  character  of  the  input  motion,  since  the  system  can  not  be  solved  for  with 
only  the  output  known. 


Tokimatsu  assumes  one  effective  soil  layer  over  bedrock.  Working  from  the  basic 
equation  linking  S-wave  velocity  and  fundamental  period  (Dobry  et  al.,  1976), 


where  H 

P 
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1 

II 
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for  the  shear  modulus  is  derived  as 

G=16p 

= thickness  of  the  soil  layer 
= mass  density  of  the  soil  layer 

[t] 

(4.14) 

= fundamental  period  of  the  soil  layer 

The  assumption  is  now  made  that  any  increase  in  the  fundamental  period  is  directly  related 
to  a decrease  in  soil  layer  stiffness: 
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where  G/G^ 


modulus  reduction  factor 

fundamental  period  of  the  soil  layer  at  lO"^  shear  strain. 


Assuming  linear  behavior  of  the  soil  mass,  a relation  for  the  shear  displacement  at 
any  depth  z from  peak  particle  velocity  and  S-wave  velocity  is  given.  This  derivation 
assumes  that  the  shear  strain  at  any  depth  is  just  a time-offset  displacement  of  the  surface 
strain.  Although  Tokimatsu  claims  now  to  be  able  to  accurately  calculate  shear  strain  at  any 
point  in  the  soil  mass  at  any  time,  he  defines  an  effective  shear  strain  to  be  85  percent  of 
the  maximum  shear  strain.  This  value  is  picked  since  for  the  first  ten  seconds  of  the 
excitation,  the  effective  strain  is  considered  to  be  85  percent  of  the  maximum  shear  strain. 

For  the  estimation  of  damping  ratio,  Tokimatsu  assumes  that  the  input  motions  from 
the  bedrock  layer  are  white  noise,  a poor  assumption  based  on  what  has  been  presented  in 
this  report.  This  assumption,  however  allows  the  authors  to  say  that  the  Fourier  spectrum 
of  the  surface  motions  has  much  the  same  shape  as  the  spectral  ratio  (i.e.  amplification  spectrum). 
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The  authors  claim  good  correspondence  between  laboratory  results  and  their 
technique.  However,  until  a more  rational  and  complete  derivation  is  given,  this  approach 
is  not  considered  useful.  There  are  too  many  theoretical  problems  and  shaky  assumptions 
involve  to  accept  it  as  presented. 

A good  example  of  the  apphcation  of  time-adaptive  modeling  of  soil  parameters  from 
earthquake  strong  motion  data  is  a paper  by  Safak  (1989).  While  the  paper  is  difficult  to 
understand,  the  application  of  advanced  system  identification  techniques  is  well  illustrated. 
In  this  paper  the  adaptive  scheme  uses  an  exponentially  decaying  "forgetting-factor"  to  take 
the  non-stationarity  into  account,  rather  than  a full  Kalman  filter  that  actually  models  the 
changing  parameters  directly. 

The  process  is  modeled  as  a single  input-single  output  system  where  the  input  and 
output  are  loiown,  but  the  auxiliary  noise  input  is  not  known  but  assumed  to  be  white.  The 
general  model  is  autoregressive-moving  average  with  a moving  average  auxiliary  input 
(ARMAX).  The  ARMA  parameters  are  estimated  by  a recursive  least  squares  scheme 
known  as  the  Recursive  Prediction  Error  Method.  Safak  shows  that  for  earthquake  ground 
motions,  a subset  of  the  ARMAX,  the  ARX  (autoregressive  with  noise)  is  most  appropriate. 
If  the  prediction  error  series  is  a stationary  Gaussian  variable  with  zero  mean,  this  method 
is  identical  to  the  Maximum  Likelihood  estimation. 

Proper  application  of  a parametric  model  to  a process  requires  determining  the 
proper  model  order.  The  simplest  criterion  is  to  use  the  order  that  minimizes  the  prediction 
error 

where  e(t)  = prediction  error  at  time  t (innovation) 

y(t)  = actual  output  at  time  t 

^(t)  = prediction  of  output  at  time  t made  at  time  t-1. 

As  discussed  in  Chapter  2,  this  estimate  can  be  checked  by  observing  whether  the  prediction 
error  time  series  is  white.  It  was  shown  that  if  all  available  information  is  retrieved  from  the 
data,  the  residuals  will  become  a white  noise  series.  However,  it  is  possible  for  the  system 
to  be  correctly  modeled  while  the  noise  is  not.  This  case  still  will  give  the  correct  transfer 
function  but  not  have  a white  residuals  series.  The  validity  of  the  model  can  still  be  checked 
by  observing  the  cross-correlation  of  the  input  and  residuals  series.  For  a valid  model , there 
will  be  no  correlation  between  the  two  since  the  model  "pulls"  all  the  information  out, 
leaving  only  the  noise.  The  final  choice  is  to  use  the  AIC,  which  maximizes  the  entropy 
between  the  model  and  actual  process. 

Safak  gives  an  example  of  the  ARX  approach  to  spectral  estimation  by  analyzing  the 
surface  strong  motion  from  the  1971  San  Fernando  earthquake.  If  the  input  white  noise 
signal  is  combined  with  the  auxiliary  white  noise  signal,  the  ARX  model  becomes  the  more 
familiar  ARMA  model.  This  method  is  identical  to  the  Burg  estimation  if  the  MA 
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parameters  (which  model  sharp  dips)  are  not  included.  As  was  shown  in  Chapter  2,  this 
method  avoids  making  the  assumption  that  the  signal  is  stationary  and  the  noise  process  is 
discrete  from  the  process  of  interest.  Figure  3.15  shows  the  resulting  instantaneous  spectral 
estimates  for  one,  five,  and  eighteen  seconds  into  the  earthquake  (note  the  log-log  scale). 
It  is  seen  that  the  fundamental  frequency  becomes  smaller  though  the  shaking. 

The  ARMAX  model  was  also  applied  to  estimating  the  parameters  of  a soil  column 
(Safak,  1989).  In  this  case  nearby  bedrock  motion  and  the  surface  motion  were  both  known 
for  a site  subjected  to  the  1985  Chilean  earthquake.  The  estimated  parameters  are 
associated  with  system  fundamental  frequency  and  damping  ratio.  In  the  frequency  domain 
the  spectrum  shows  frequency-dependant  amplification.  The  method  used  allows  the  process 
parameters  to  be  estimated  throughout  time  so  that  the  change  in  the  system  can  be 
monitored.  In  this  paper  only  the  results  for  eighteen  seconds  into  the  temblor,  the  coda, 
is  reported.  The  results  for  the  first  mode  are  a natural  frequency  of  1.35  Hz  and  a damping 
ratio  of  11.9  percent. 

The  results  in  this  paper  are  arrived  at  through  a deliberate,  rational  approach  rather 
than  trial-and-error,  ad  hoc  methods.  Obvious  areas  of  improvement  would  be  to 
incorporate  some  of  the  techniques  developed  by  Gersch  and  his  cohorts.  The  author 
himself  points  out  that  the  ARMAX  model  did  not  model  the  noise  very  well,  and  suggests 
the  use  of  the  Box-Jenkins  model  where  the  noise  is  modeled  with  it’s  own  ARMA 
polynomial.  Safak  also  points  out  that  the  application  of  a Kalman  filter  would  take  the 
non-stationarity  more  rationally  into  account. 
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Amplitude  (cm/s) 


Fig.  3.15  Amplitude  of  the  transfer  functions  of  the  1971  San  Fernando  temblor  modeled 
as  an  ARMA  (8,7)  process  at  times  t = 1,  5,  and  18  seconds  (from  Safak,  1989). 
Note  log-log  scale. 
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CHAPTER  4 CONCLUSIONS  AND  SUMMARY 


4.1  A Framework  of  Knowledge 

4.1.1  Estimation  Techniques 

For  the  system  identification  problems  of  interest  to  this  report,  the  parameters  to 
be  found  are  system  damping  and  resonant  frequency.  The  resonant  frequency  is  associated 
with  the  peak  of  the  amplitude  spectrum.  It  has  been  seen  that  a single  peak  is  relatively 
easy  to  identify,  especially  if  the  proper  model  has  been  used  for  the  process.  For  the  very 
large  resonance  peak  of  the  principle  harmonic,  virtually  any  method  will  give  a close  result, 
even  for  a relatively  non-stationary  case.  However,  every  method  distorts  the  shape  of  the 
peaks,  the  usual  source  of  the  damping  estimate.  It  is  therefore  to  be  expected  that  system 
identification  will  give  good  estimates  of  modal  frequencies,  but  poor  estimates  of  damping. 

A study  was  done  by  Gersch  (1974)  in  order  to  determine  the  greatest  degree  of 
accuracy  with  which  a proper  order  ARMA  model  can  estimate  the  damping  ratio  and 
natural  frequency  of  a structure,  using  the  Maximum  Likelihood  method.  He  notes  that  as 
the  number  of  data  points  (N)  becomes  large,  the  estimate  approaches  the  actual  model  and 
approaches  the  Cramer-Rao  lower  bound  of  variance.  The  accuracy  is  given  as  the 
coefficient  of  variation  — standard  deviation  divided  by  the  mean.  For  both  parameters,  the 
coefficient  of  variation  is  inversely  proportional  to  JN  and  length  of  the  period  sampled,  and 
relatively  insensitive  to  noise.  For  one  thousand  data  points,  the  coefficient  of  variation  was 
less  than  0.01  for  natural  frequency,  but  greater  than  0.2  for  damping  ratio;  this  is  an  order 
of  magnitude  difference. 

As  a note  of  caution,  while  the  more  complicated  models  encompass  more  of  a 
process  in  a rational  manner,  the  values  needed  as  input  must  often  be  assumed,  making  the 
results  less  meaningful.  The  results  of  the  more  complicated  models  should  be  compared  to 
those  of  simpler,  more  intuitive  calculations  that  allow  the  engineer  to  maintain  a "feel"  for 
the  validity  of  the  results  throughout  the  evaluation  process.  In  the  words  of  G.  E.  P.  Box, 
"All  models  are  incorrect,  but  some  are  more  useful  than  others." 


4.1.2  Conclusions  and  Recommendations 

Theoretical  and  practical  considerations  have  shown  that  transforming  a data  series 
into  the  frequency  domain  is  not  a straight-forward  process.  The  results  of  a simple  DFT 
has  variance  equal  to  the  mean,  and  is  very  severely  corrupted  by  leakage.  This  has  resulted 
in  the  tremendous  amount  of  ad  hoc  attempts  to  improve  on  the  bias  and  uncertainty 
problems,  lumped  in  this  paper  as  the  "classical  approach".  This  presentation  has  not  gone 
into  detail  about  these  methods,  except  to  become  aware  of  the  inherent  problems.  A 
thorough  summary  and  comparison  of  almost  all  these  non-parametric  methods,  and 
parametric  methods,  can  be  found  in  Kay  and  Marple  (1981).  Included  Kay  and  Marple’s 


paper  is  a shocking  comparison  of  the  results  of  all  the  methods,  each  giving  results 
unrelated  to  the  others,  and  all  different  from  the  theoretical  frequency  spectrum. 

These  problems  were  directly  addressed  by  Thomson  (1982)  in  deriving  the  multi- 
taper method.  It  is  believed  that  this  approach  is  the  ’'best"  non-parametric  spectral 
approach  if  proper  pre-treatment  of  the  signal  is  done.  These  same  considerations  are 
addressed  from  a philosophical  point  of  view  by  many  papers  of  John  Tukey  (e.g.  Tukey, 
1984;  Brillinger  and  Tukey,  1984).  Blindly  applying  the  FFT  to  one’s  data  is  a very 
dangerous  thing  to  do,  and  the  cautions  and  insights  of  Tukey  should  be  taken  to  heart 
before  starting. 

If  the  process  being  studied  can  be  modeled  as  an  AR  or  ARMA  process,  then  the 
parametric  approach  is  the  best  method  to  characterize  the  system.  The  Burg  method  is 
ideal  for  short,  relatively  stationary  data.  The  growing  family  of  adaptive  and  Kalman  filters 
are  proving  themselves  with  non-stationary  processes.  The  AR  model  was  shown  to  be 
initially  derived  for  a SDOF  oscillator,  and  the  ARMA  model  can  be  derived  directly  from 
the  differential  equation  of  motion  for  an  N-degree-of-freedom  system,  with  the  damping 
ratio  and  resonant  frequency  the  model  parameters  (e.g.  Gersch  and  Luo,  1970).  A 2n-2n 
ARMA  model  is  therefore  a valid  model  for  a layered  soil  system,  or  soil-structure 
interaction  problem. 

The  examination  of  the  various  attempts  to  characterize  soil  properties  through 
analysis  of  the  response  to  earthquake  excitation  has  illustrated  some  of  the  main  pitfalls  and 
advantages  of  using  the  system  identification  approach.  One  important  problem  in 
evaluating  the  various  methods  is  that  there  is  no  "correct"  value  against  which  to  compare 
the  results.  The  mechanical  engineer  has  the  advantage  of  being  able  to  construct  a system 
similar  to  that  being  tested,  with  known  parameters  against  which  to  test  the  method.  The 
geotechnical  engineer  never  has  this  luxury.  This  limitation  highlights  one  of  the  strengths 
of  the  SI  approach,  since  there  is  no  other  way  to  actually  measure  the  in  situ  fundamental 
frequency  and  damping  ratio  for  strains  even  remotely  approaching  those  encountered 
during  earthquake  loading. 

The  geotechnical  community  has  not  utilized  the  methods,  approaches,  or  warnings 
discussed  by  the  system  identification  community  and  presented  in  Chapter  2.  The  methods 
used  to  date  make  no  acknowledgement  to  the  non-stationarity  of  the  signals,  or  that  the 
Fourier  spectral  estimate  may  have  any  limitations  or  uncertainties.  The  exceptions  to  this 
are  the  researchers  using  forced  vibration  and  fitting  impedance  functions. 

The  forced  vibration-impedance  function  method  is  very  attractive  since  the 
investigators  have  control  over  the  input  signal.  Besides  simplifying  the  calculations,  this 
avoids  having  to  wait  many  years  for  an  earthquake  to  occur  at  a given  site.  The  approach 
is  also  straight-forward  in  concept.  The  major  problem  seems  to  be  an  inability  of 
generating  enough  energy  to  involve  a deep  column  of  soil.  The  strains  induced  in  the  soil 
will  be  in  the  low  range  for  the  same  reason. 
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It  is  recommended  that  the  complete  input-output  data  sets  from  the  Wildlife  site  and 
Lotung  be  analyzed  using  the  tools  provided  by  the  researchers  in  control.  The  improved 
segmentation  scheme  of  Gersch  and  Brotherton  (1982)  is  very  appealing  if  the  process 
parameters  are  not  believed  to  change  during  excitation.  For  situations  where  the  soil  is 
believed  to  have  undergone  degradation,  methods  taken  after  ICitagawa  and  Gersch  (1985) 
seems  to  offer  the  most  promise.  In  any  event,  application  of  AR  and  ARMA  spectral 
estimators  will  yield  spectral  estimates  with  higher  certainty  than  Fourier  analysis. 
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